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ABSTRACT 


We  develop  a  numerically  stable  algorithm  for  electromagnetic  wave  propagation  through 
planar  stratified  media.  This  algorithm  is  implemented  in  a  modern  programming  language 
and  is  suitable  for  the  study  of  such  applications  as  Anderson  localization  and  perfect  lens- 
ing.  Our  algorithm  remains  numerically  stable  even  in  the  presence  of  large  absorption. 
Furthermore,  in  the  context  of  the  linear  response  laws  and  causality,  we  analyze  a  vanish¬ 
ing  absorption  approximation,  which  is  commonly  used  in  wave  scattering  problems.  We 
show  that  it  is  easy  to  violate  causality  in  the  frequency-domain  by  making  the  vanishing 
absorption  approximation. 

We  also  develop  an  orders-of-scattering  approximation,  termed  “screened  cylindrical 
void/core”  (SCV)  approximation,  for  wave  scattering  from  a  large  host  cylinder  contain¬ 
ing  N  eccentrically  embedded  core  cylinders.  The  SCV  approximation  is  developed  via 
separation  of  variables  and  a  cluster  T-matrix.  We  establish  the  limitations  of  the  SCV 
approximation  and  it  is  in  good  agreement  with  the  numerically-exact  solution.  Fur¬ 
thermore,  we  illustrate  that  the  large  host  cylinder  model  with  N  cylindrical  inclusions  can 
be  used  to  theoretically  and  experimentally  investigate  strong  multiple  scattering  effects  in 
random  media,  such  as  Anderson  localization. 
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CHAPTER  1 


INTRODUCTION 

It  has  over  150  years  since  James  Clerk  Maxwell  published  the  much-celebrated  equations 
that  now  bear  his  name.  He  probably  could  not  have  imagined  that,  in  the  21st  century,  his 
equations  would  be  utilized  to  study  various  wave  scattering  phenomena.  The  main  goal  of 
this  dissertation  is  to  use  Maxwell’s  equations  to  design  multiple  scattering  experiments  that 
can  be  realized  in  the  Mesoscopic  Physics  Laboratory  (MPL)  at  the  Colorado  School  of  Mines. 
These  modeled  experiments  not  only  need  to  be  theoretically  sound,  but  they  also  need  to 
be  computationally  and  experimentally  verihable.  At  the  MPL,  we  have  the  rare  ability 
to  measure  both  the  amplitude  and  the  phase  of  electric  helds  at  frequencies  approaching 
1  THz  with  the  AB  Millimetre  millimeter /sub-millimeter  Vector  Network  Analyzer  (VNA), 
see  Figure  1.1.  This  unique  ability  allows  us  to  study  one  of  the  most  fascinating  multiple 
scattering  phenomena,  namely,  the  localization  of  electromagnetic  waves.  Fifty  years  after 
the  publication  of  Anderson’s  seminal  work  [1],  localization  continues  to  be  a  thriving  area 
of  research  in  theoretical  and  experimental  physics  [2-7]. 

Localization,  a  coherent  multiple  scattering  effect,  has  been  observed  in  a  variety  of 
classical  and  quantum  mechanical  systems.  To  gain  a  pictorial  understanding  of  how  a 
wave  may  become  spatially  localized,  consider  a  wave,  incident  on  a  medium 

composed  of  N  identical  and  lossless  scatterers  randomly  distributed  in  space.  The  total 
wave,  is  the  sum  of  the  incident  and  partially  scattered  waves,  i.e.,  (r)  = 

_l_  ^^^7/W(r),  where  U^^\r)  is  the  wave  scattered  by  the  Rh  scatterer.  These 
waves  may  interfere  constructively  in  some  region  of  space  (see  Figure  1.2)  and  destructively 
in  other  regions  of  space,  thus  localizing  in  space. 
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Figure  1.1:  John  A.  Scales  (top)  and  Philippe  Goy  (bottom)  calibrating  MPL’s  AB  Millime¬ 
tre  VNA,  which  Goy  designed  and  build  [8]. 


Figure  1.2:  A  schematic  representation  of  constructive  interference  of  three  partially  scat¬ 
tered  waves  is  shown.  The  scatterers  are  denoted  by  disks  and  the  region  of  constructive 
interference  is  highlighted. 

1.1  One-dimensional  scattering  and  absorption  free  media 

The  hrst  use  of  the  VNA  at  the  MPL  was  to  study  localization  in  planar  stratihed 
media  (see  Figure  1.3)  composed  of  100  randomly  shuffled  quartz  and  Teflon  layers  [7].  In 
the  course  of  this  study,  it  became  apparent  that  the  standard  transfer  matrix  algorithm 
[9-11]  used  to  compute  the  electromagnetic  field  inside  the  multilayer  stack  is  numerically 
unstable.  The  root  cause  of  the  numerical  instability  is  that  the  transfer  matrix  contains 
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exponentially  increasing  and  decreasing  terms  when  absorption  is  non  zero.  In  Chapter  2,  we 
present  a  nnmerically  stable  algorithm  for  a  multilayer  stack,  which  we  have  developed  and 
implemented  in  a  modern  object-oriented  programming  language.  Furthermore,  we  discuss 
some  pathological  cases  that  arise  in  the  limit  of  vanishing  absorption.  These  cases  are 
important  as  they  vividly  demonstrate  that  there  is  no  physical  principle  that  can  be  used 
to  differentiate  between  a  right-handed  material  and  a  left-handed  material  if  absorption  is 
assumed  to  be  zero. 


incident  I 
wave\  I  ^ 


\ 

I 

Figure  1.3:  The  cross-sectional  view  of  the  planar  stratified  media  (multilayer  stack)  is  shown. 


The  common  zero  absorption  approximation  is  further  investigated  in  Chapter  3  from  a 
more  general  point  of  view.  In  particular,  a  connection  between  absorption  and  causality 
(the  effect  cannot  precede  the  cause)  is  established  in  the  context  of  linear  response  laws. 
Furthermore,  through  the  use  of  the  formal  theory  of  tempered  distributions,  we  estab¬ 
lish  significant  differences  between  the  frequency-domain  and  monochromatic  time-domain 
Maxwell  equations.  This  is  important  because  although  all  the  scattering  problems  in  this 
dissertation  are  formulated  in  the  frequency-domain,  the  actual  VNA  measurements  are  per¬ 
formed  in  the  time-domain  with  a  monochromatic  incident  field.  As  shown  in  Chapter  3, 
confusing  the  two  domains  leads  to  an  apparent  violation  of  causality. 

^It  is  customary,  but  historically  inaccurate,  to  attribute  the  first  theoretical  consideration  of  the  left-handed 
material  to  Veselago  [12].  In  fact,  as  Agranovich  and  Gartstein  point  out,  such  a  material  was  theoretically 
considered  much  earlier  by  Mandel’shtam,  see  [13,  §2.1]  and  references  therein. 
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1.2  Large  host  cylinder  with  holes 


At  the  MPL,  we  are  currently  fabricating  a  model  of  a  millimeter-wave  random  medium 
from  a  large  cylinder  of  Teflon  (ultra  low-loss  material)  with  approximately  5000  holes  drilled 
in  a  random  pattern  using  a  table-top  three- axis  CNC  milling  machine.  The  cylinder  is 
110  mm  in  length  and  has  a  diameter  of  152.4  mm.  The  holes  are  0.505  mm  in  diameter 
and  are  separated  by  a  distance  of  0.101mm,  giving  us  roughly  two/three  scatterers  per 
wavelength,  as  required  by  the  loffe-Regel  criterion  [14,  §7.4.4]  for  localization.  The  holes 
are  large  enough  to  allow  for  the  possibility  of  placing  an  intensity  detector  in  them,  thereby 
measuring  inside  the  sample.  Our  model,  shown  in  Figure  1.4,  can  be  illuminated 

from  the  side  to  make  a  two-dimensional  system  or  from  the  end  to  study  transverse  localiza¬ 
tion.  Furthermore,  by  illuminating  the  sample  from  the  side  and  putting  it  on  a  rotational 
stage,  we  can  generate  essentially  arbitrary  realizations  of  the  same  random  disorder.  In  this 
dissertation,  we  will  only  consider  illumination  from  the  side;  however,  it  is  important  to 
realize  that  the  above  model  is  more  versatile  than  considered  here. 


Figure  1.4:  Artistic  rendering  (not  to  scale)  of  our  millimeter- wave  random  medium  model. 
The  holes  are  shown  in  a  simple  pattern  for  illustration  purposes  only. 
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In  order  to  ensure  that  the  VNA  has  enough  sensitivity  to  detect  the  presence  of  such 
small  holes,  we  measured  the  total  held  in  the  forward  direction  as  a  function  of  distance 
with  and  without  one  concentric  hole.  The  results  of  this  experiment  are  shown  in  Figure  1.5, 
from  which  we  clearly  see  that  the  VNA  can  unambiguously  detect  the  presence  of  a  small 
hole  inside  a  large  Tehon  cylinder. 


Figure  1.5:  The  measured  amplitude  and  phase  of  the  total  transmitted  held  at  160  GHz  are 
shown. 

We  believe  that,  in  order  to  learn  the  most  from  the  aforementioned  experiments,  it  is 
necessary  to  have  a  sound  theoretical  model  for  them.  The  theoretical  model  should  be 
realistic  but  computationally  tractable,  without  many  crude  numerical  approximations.  It 
is  possible  to  model  our  ensemble  of  scatterers  using  integral  equation  methods  [15-17]  or 
coupled-dipole  methods  [18,  19],  but  such  methods  are  not  computationally  feasible  because 
they  yield  a  linear  system  of  equations  that  is  too  large.  For  example,  it  is  common  practice 
to  use  supercomputers  with  the  coupled-dipole  method  to  model  even  a  few  scatterers  of 
wavelength  size  [19].  Instead,  in  Chapters  4-6,  we  have  developed  an  orders-of-scattering 
approximation  based  on  separation  of  variables  and  a  cluster  T-matrix  [20,  21]  [22,  §5.9]  [23, 
§7.11].  This  development  is  done  in  stages,  as  described  below. 
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In  Chapter  4,  we  develop  an  approximation,  termed  the  screened  cylindrical  void/core 
(SCV)  approximation,  for  a  large  cylinder  with  one  concentric  hole.  The  SCV  approximation 
is  hrst  derived  in  a  physically  intnitive  manner,  and  then  in  a  mathematically  rigorous  man¬ 
ner.  In  Chapter  5,  the  SCV  approximation  is  extended  to  the  case  of  a  non-concentric  hole 
and  we  explicitly  show  that  the  SCV  approximation  can  be  viewed  as  an  orders-of-scattering 
approximation.  Furthermore,  the  limitations  of  the  SCV  approximation  are  discussed  and 
demonstrated  with  explicit  numerical  examples.  Finally,  in  Chapter  6,  we  generalize  the  SCV 
approximation  to  N  non-concentric  holes  inside  the  large  cylinder  via  the  cluster  T-matrix 
and  suggest  avenues  for  future  research. 
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2.1  Abstract 

The  A-matrix  algorithm  for  the  propagation  of  an  electromagnetic  wave  through  planar 
stratihed  media  has  been  implemented  in  a  modern  object-oriented  programing  language. 
This  implementation  is  suitable  for  the  study  of  such  applications  as  the  Anderson  localiza¬ 
tion  of  light  and  super- resolution  (perfect  lensing).  For  our  open-source  code  to  be  as  useful 
as  possible  to  the  scientihc  community,  we  paid  particular  attention  to  the  pathological  cases 
that  arise  in  the  limit  of  vanishing  absorption. 

2.2  Introduction 

Electromagnetic  wave  propagation  through  planar  stratihed  media  (multilayer  stack)  is  a 
century  old  problem  in  physics  [1,  2].  It  may  be  somewhat  surprising  that  it  is  still  relevant 
today.  In  fact,  it  has  only  relatively  recently  been  discovered  that  the  transmission  and 
rehection  coefficients  for  a  multilayer  stack  may  be  written  down  without  any  computations 
by  using  a  complex  version  of  the  elementary  symmetric  functions  [3,  4].  It  has  also  been 
recently  discovered  that  the  complex  rehection  coefficients  follow  the  generalized  version  of 
the  composition  law  used  to  add  parallel  velocities  in  the  theory  of  special  relativity,  see 
[5,  6]  and  Refs,  within.  It  is  possible  to  use  the  aforementioned  properties  to  formulate  a 
numerical  wave  propagation  algorithm  in  planar  stratihed  media  as  was  done  in  [7],  yet  the 


resulting  algorithm  appears  to  be  numerically  unstable.  The  more  traditional  approach  of 
the  late  1940s,  namely,  the  transfer  matrix  algorithm  [8-11],  is  also  numerically  unstable. 
Both  algorithms  are  numerically  unstable  because  they  contain  exponentially  increasing  and 
decreasing  terms,  see  Section  2.6.  There  also  exists  an  i?-matrix  algorithm  [12-15],  but  it 
is  only  conditionally  stable  (for  reasons  different  from  above)  [12,  15].  We  use  a  simple 
version  of  the  S'-matrix  algorithm,  which  is  numerically  stable  [15-19].  Before  considering 
the  details  of  the  S'-matrix  algorithm  and  the  need  for  its  open-source  implementation  in 
a  modern  object-oriented  language,  we  briefly  mention  some  of  the  current  applications  we 
had  in  mind  when  we  wrote  the  code. 

In  1968,  Veselago  [20]  considered  a  hypothetical  non-active  material  in  which  the  real 
parts  of  the  permittivity  and  permeability  are  simultaneously  negative;  we  refer  to  such 
a  material  as  a  left-handed  material  (LHM),  but  it  is  also  known  as  a  negative  refractive 
material.  It  was  only  in  the  early  2000’s  that  such  an  artihcial  material  was  fabricated  [21,  22], 
leading  to  an  explosion  of  papers  on  the  LHM,  see  [23]  and  Refs,  within.  One  of  the  intriguing 
properties  of  the  LHM  is  the  ability  to  image  with  a  sub-wavelength  image  resolution  (super¬ 
resolution  if  you  will),  which  has  been  proposed  and  studied  in  the  context  of  a  multilayer 
stack  [24,  25].  Another  general  area  of  application  is  the  Anderson  localization  of  light 
[26,  27],  which  has  been  studied  both  theoretically  and  experimentally  by  Scales  et  ah  [28], 
who  considered  wave  propagation  at  normal  incidence  through  a  multilayer  stack  made  of 
quartz  and  Teflon  wafers.  The  effects  of  total  internal  reflection  on  light  localization  in  a 
random  multilayer  stack  at  oblique  incidence  have  also  been  studied  under  the  assumption 
of  complete  phase  randomization  [29]  along  with  the  effects  of  the  LHM  on  localization  [30] . 
Other  applications  include  the  study  of  asymmetrical  properties  of  light  in  a  Fabry-Perot 
interferometer  [31,  32]. 

In  all  of  the  above  applications,  the  S'-matrix  algorithm  was  or  could  have  been  used; 
however,  to  the  best  of  our  knowledge,  an  open-source  and  object-oriented  implementation 
of  the  S'-matrix  algorithm  suitable  for  the  LHM  as  well  as  the  right-handed  material  (RHM) 
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(where  the  real  parts  of  the  permittivity  and  permeability  are  not  simultaneously  negative) 
is  currently  unavailable.  Almost  certainly,  there  are  many  “in-house”  implementations  of 
some  version  of  the  algorithms  discussed  above  being  passed  around  among  colleagues.  We 
suspect  that  some  users  of  these  “in-house”  algorithms  may  be  unaware  of  the  numerical 
stability  issues  and  of  pathological  cases  where  the  numerical  implementation  is  not  clear, 
as  discussed  in  Section  2.4.  Moreover,  in  the  context  of  reproducibility  of  scientihc  work,  it 
is  important  to  have  an  open-source  and  publicly  available  implementation. 

This  paper  is  self-contained  as  much  as  possible  in  order  for  our  implementation  of  the 
A-matrix  algorithm  to  be  useful  to  the  widest  possible  scientihc  community.  We  also  point 
out  the  benehts  and  drawbacks  of  using  a  high-level  programing  language  called  Python  for 
implementing  our  code,  see  Section  2.10. 

2.3  Background 

The  source-free  macroscopic  Maxwell  equations  with  assumed  harmonic  time  dependence, 
exp  (— icat),  in  the  Systeme  International  (SI)  unit  system,  at  every  ordinary  point  in  space, 
are: 


V-D  =  0,  V-B  =  0,  (2.1a) 

V  X  E  =  icaB,  V  x  H  =  — icaD,  (2.1b) 

where  E  is  the  electric  held,  D  is  the  displacement  held,  B  is  the  magnetic  held,  H  is 
the  magnetic  intensity,  and  u  is  the  angular  frequency.  By  an  ordinary  point  in  space,  we 
mean  a  point  in  space  in  whose  “neighborhood”  the  physical  properties  of  the  medium  are 
continuous.  Thus,  strictly  speaking,  one  cannot  apply  Maxwell’s  equations  at  a  surface  that 
separates  two  physically  diherent  media.  If  the  medium  is  isotropic  and  homogeneous,  then 
D  =  eE  and  B  =  //H,  where  e  and  ^  are  the  permittivity  and  the  permeability,  respectively. 
Permittivity  must  satisfy  the  Kramers-Kronig  relations  and  is  therefore  a  complex-valued 
function  of  angular  frequency.  The  same  is  true  for  permeability.  Thus,  in  general,  we  have 
e  =  e(a;)  G  C  and  fi  =  /i(a;)  G  C. 
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The  source-free  macroscopic  Maxwell  equations  are  first-order  linear  partial  differential 
equations  (PDEs)  that  must  be  supplemented  by  some  boundary  conditions.  The  conven¬ 
tional  boundary  conditions  for  a  source-free  interface  separating  two  media  (1  and  2)  are: 


n- (D(2) =  0,  n- (B(2) -BW)  =0,  (2.2a) 

n  X  (E(2)  -  =0,  n  X  (H^^)  -  =  0,  (2.2b) 

where  n  is  a  unit  normal  to  the  interface,  and  the  superscript  on  the  helds  indicates  from 
which  medium  the  interface  is  approached. 

Taking  the  curl  of  (2.1b),  then  simplifying  the  result  using  the  V  x  (V  X  A)  =  V(V-A)- 
V^A  vector  identity  and  (2.1a),  we  obtain  the  vector  Helmholtz  equation  within  each  layer 

(V'  +  A:')|h}=0>  (2-3) 

where  k  is  the  complex  wavenumber,  and  In  general,  k‘^  ^  kk*,  where  *  denotes 

the  complex  conjugate,  and  the  computation  of  k  from  k'^  must  be  done  with  extreme  care. 
For  example,  the  permittivity  and  permeability  for  an  absorbing  material  are  taken  to  be 
e  =  e'  -|-  ie"  and  jj,  =  jj,'  +  i/i",  respectively,  where  {e',  p'}  G  M,  {e",  p"}  G  M"*"  and  M"*"  denotes 
the  positive  real  numbers.^  Let  e  =  |e|e‘®'  and  fi  =  |/i|e‘^^,  where  {9^,9^}  G  [0,7r].^  Then 

k  =  efxoo 

. -  ■(  “hSTTTi  N 

k  =  A/|e||/i|a;e  V  ^  n  =  0, 1,  (2.4) 

where  a;  >  0.  The  choice  of  the  root  in  (2.4)  is  dictated  by  the  physical  requirement 
that,  in  an  absorbing  medium,  the  wave  must  decay  and  not  exponentially  grow.  Let  k  = 
k'  +  ik",  {k',  k”}  G  M.  Without  loss  of  generality,  consider  a  plane  wave  propagating  in  the 
positive  x-direction;  then,  we  have  ^  Therefore,  k”  must  be  greater 

^For  the  exp  (-|-iwt)  time  dependence,  e  =  e'  —  ie",  g  =  g'  —  ig",  where  {e',  g'}  €  M,  fe",  g"}  G  M"*". 

^We  always  mean  the  positive  square  root  of  x  when  we  write  \/x,  where  x  G  K“'".  The  fundamental  issue 
with  the  w  =  mapping  is  that  the  “square  root”  function  has  branch  points  at  z  =  0  and  z  =  oo  and 
thus  must  have  a  branch  cut  connecting  the  two  branch  points,  see  [33,  vol.  1,  Sec.  54]. 
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than  zero  in  order  for  the  wave  to  decay  in  the  positive  x-direction. 

2.3.1  Pathological  cases  at  normal  incidence 

In  the  case  of  a  perfect  dielectric  (e"  =  0  and  jj,"  =  0),  the  rule  for  choosing  a  physically 
appropriate  root  in  (2.4)  may  be  established  by  taking  the  limit  as  absorption  goes  to  zero. 

Consider  an  almost  perfect  dielectric  made  of  the  RHM.  Let  e  =  |e|e'®%  p  =  |p|e‘^^,  where 
9e  and  9^  are  inhnitesimally  small  positive  numbers,  then  "  ^  ^  tt  and  -^y-^+7r  >  tt.  Thus, 
we  must  choose  the  n  =  0  root  in  (2.4),  i.e.,  k  =  A/|e||/i|e  v  ^  'u.  In  the  case  of  a  truly 

perfect  dielectric  (at  hxed  frequency),  we  may  take  the  limit  as  9^  and  9^  approach  zero  to 
obtain  k  =  \/\e\ \p,\u. 

In  the  case  of  an  almost  perfect  dielectric  made  of  a  LHM:  Let  e  =  |e|e‘®%  /i  =  |/i|e‘®'", 

0  -\-0  0  I  0 

where  9^  and  9^  are  slightly  less  than  tt,  then  "  ^  ^  <  tt  and  +  tt  >  tt.  Thus,  we 

-  j/ y+y  N 

must  again  choose  the  n  =  0  root  in  (2.4),  i.e.,  k  =  A/|e||/i|e  v  ^  'u.  For  a  truly  perfect 

dielectric  (at  fixed  frequency),  we  may  take  the  limit  as  9^  and  9^  approach  tt  to  obtain 
k  =  Vkllhl  =  —  A/|e| \ia\uj.  Notice  that  for  the  LHM  with  zero  absorption,  k  <  0,  and 
for  the  RHM  with  zero  absorption,  k  >  0. 

2.4  Wave  propagation  in  stratified  media 


Consider  the  three-dimensional  space  divided  into  p  -|-  1  regions.  The  regions  are  inhnite 
in  the  yz-pAane,  see  Figure  2.1. 

The  interfaces  separating  the  regions  are  assumed  to  be  perfectly  planar  [yz-plane).  The 
regions  ^  =  0,  ...,p  —  1  are  assumed  to  be  isotropic  and  homogeneous  with  a  complex 
permittivity,  e^,  and  complex  permeability,  The  region  p  is  assumed  to  be  isotropic  and 
homogeneous  with  real  permittivity,  Cp,  and  real  permeability,  Pp.  In  other  words,  we  have 
{e^,  ye}  G  C  for  £  =  0, . . . ,  p  —  1  and  {e^,  Pp]  G  M. 

A  monochromatic  plane  wave  in  the  Hh  region  is  given  by 


Ee{r,t) 

He{r,t) 


Ee 


—  J  I  pi(krr-a;t) 


=  0, 


,p, 


(2.5) 
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Figure  2.1:  The  cross-sectional  view  of  the  multilayer  stack  is  shown.  The  multilayer  stack 
consists  of  p  -|-  1  regions  made  of  a  RHM.  A  parallel  polarized  wave  is  incident  from  a  semi- 
inhnite  ambient  medium  (region  p).  The  origin  of  the  coordinate  system  is  set  on  the  planar 
interface  separating  regions  p  and  p  —  1.  The  0th  region  is  a  semi-inhnite  substrate. 


where  r  =  xH  +  yy  +  zz,  {E^,  are  the  complex  vector  amplitudes,  x  -|-  ky^i  y  -|- 

kz/z  is  the  complex  wavevector.  It  is  clear  that  (2.5)  satishes  (2.3)  if 

=  kle  +  kl^  +  klg  =  kj  =  (2.6) 

Without  loss  of  generality,  we  can  set  kz,e  =  0  because  we  can  always  rotate  the  coordinate 
system  so  that  the  p-axis  is  parallel  to  the  part  of  the  k  vector  that  lies  in  the  pz-plane, 
see  Figure  2.1.^  The  solution  given  by  (2.5)  in  each  region  must  also  satisfy  the  boundary 
conditions  given  by  (2.2).  Substituting  (2.5)  into  (2.2)  yields, 

ky,p  =  ky^e,  i  =  0, . . .  ,p  -  1,  (2.7) 

where  ky^p  G  M  because  we  have  assumed  that  the  region  p  has  real  permittivity  and  per¬ 
meability.  Therefore,  from  (2.7)  we  have  ky^i^  G  M,  but  note  that  in  general,  k^/  G  C  for 


^We  could  have  chosen  to  set  =  0,  and  then  rotated  the  coordinate  system  so  that  the  2;-axis  is  parallel 
to  the  part  of  the  k  vector  that  lies  in  the  yz-plane.  The  point  is  that  k  can  always  be  made  into  a 
two-dimensional  vector. 
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i  =  0, . . .  ,p  —  1.  Using  (2.6)  and  (2.7)  yields 

kx,e  =  with  >  0  (2.8) 

for  £  =  0, . . .  ,p,  where  Im  denotes  the  imaginary  part,  and  the  root  choice,  >  0,  is 

dictated  by  the  decaying  wave  requirement,  see  Section  2.3. 

2.4.1  Pathological  cases  at  oblique  incidence 

It  is  clear  from  (2.8)  that  if  e"  =  0,  /i"  =  0  and  >  ky  ^,  then  the  root  choice  is  not 

resolved  by  the  \m[kx^(\  >  0  requirement.  In  order  to  resolve  the  root  choice,  we  proceed  by 
taking  a  limit  as  absorption  goes  to  zero  just  as  we  did  in  Section  2.3.1.  For  the  RHM,  let 

=  |/i£|e‘^'"«  and  for  the  LHM,  let  /i^  =  where  9^^ 

and  6*^^  are  inhnitesimally  small  positive  numbers.  Then  ^  can  be  approximately  written 
as  kl^^  ~  |R|  where  0  <  7  tt,  hm|^„  ^//|_j^g7  =  0,  and  the  positive  (negative)  sign  in 
the  exponential  corresponds  to  the  RHM  (LHM).  Thus,  we  have 

Im  [kx,e\  =  a/I^  I  sin 

where  it  is  clear  that  for  the  RHM  (LHM)  the  hrst  (second)  root  must  be  chosen  in  order 
for  Im  [kx/]  >  0.  Therefore,  if  e"  =  0,  /i"  =  0  and  eiPiu"^  >  ky  p,  then  for  the  RHM  we  have 

kx,e  =  Ihd  -  k^^p,  and  for  the  LHM  we  have  kx,e  =  -\J\ei\  \pn\ -  k^p. 

2.4.2  Origin  and  numerical  treatment  of  the  pathologies 

The  limiting  procedure  carried  out  in  Section  2.3.1  and  2.4.1  appears  to  be  reasonable, 
but  unfortunately,  it  is  also  not  physically  attainable,  even  in  principle!  If  we  view  e(a;)  and 
/i(a;),  where  u  =  u'  +  ica",  in  the  context  of  the  Kramers-Kronig  relations,  then  e(cj)  and 
/u(cv)  are  analytic  functions  in  the  upper-half  ca-plane.  Furthermore,  it  can  be  shown  that 
€(lv)  and  ju(cv)  are  never  purely  real  for  any  hnite  lv  except  for  lv'  =  0  (positive  imaginary 
axis),  e.g.,  see  [34,  Section  123]  and  [35,  Section  82].  Therefore,  the  common  practice  of 
replacing  e'  +  ie"  by  e'  and  pt!  +  ip"  by  p'  even  in  an  inhnitesimally  small  u'  interval  cannot 
be  justihed.  Moreover,  by  considering  the  global  behavior  of  kx,e  it  can  be  shown  that  for  a 
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non-active  medium  is  never  zero  [36].  However,  we  see  from  (2.8)  that  kx/,  for  any  ^  ^  p 
may  be  equal  to  zero  if  and  /i£  are  purely  real.  Of  course,  this  case  only  occurs  when  the 
angle  of  incidence  precisely  equals  one  of  the  critical  angles,  and  from  the  global  properties 
of  e  and  p  we  see  that  such  angles  cannot  exist. 

The  above  discussion  suggests  that  the  pathological  cases  only  occur  in  an  unphysical 
approximation,  i.e.,  e  ~  e'  and  p  ~  p' .  In  our  numerical  code,  the  user  may  select  how  to 
deal  with  the  pathologies  from  the  following  two  schemes: 

1.  If  a  region  contains  purely  real  permittivity  and  permeability,  then  the  real  permittivity 

and  permeability  are  replaced  by  a  slightly  absorbing  permittivity  and  permeability, 
respectively,  i.e.,  for  £  7^  p,  — )■  -|-  ie"  and  p'^  ^  p'^  +  ip”,  where  e”  and  p”  are  small 

positive  numbers. 

2.  If  a  region  contains  purely  real  permittivity  and  permeability,  then  the  k^/  is  computed 
as  describe  in  Section  2.3.1  and  2.4.1.  If  this  scheme  is  chosen,  then  the  code  may 
produce  erroneous  results  at  or  very  near  the  critical  angles. 

2.5  Polarization 


The  most  general  polarization  state  is  an  elliptical  polarization  state.  However,  there 
is  no  need  to  consider  this  general  case  because  an  elliptical  polarization  state  can  always 
be  decomposed  into  a  linear  combination  of  two  linearly  independent  polarization  states, 
namely,  the  parallel  polarization  state  and  the  perpendicular  polarization  state.  In  what 
follows,  it  is  convenient  to  express  E£(r,  t)  and  H£(r,  t)  in  terms  of  each  other  by  substituting 
(2.5)  into  (2.1b)  (with  Dg  =  e^E^)  and  using  the  vector  identity 


V  X 


Ei{r,t) 

H,(r,t) 


=  ik/i  X 


Eg 


to  obtain 


15 


Ei{r,t) 


ke  X  H£(r,t) 


(2.9a) 


eeoj 

kj  X  Ei{r,t) 
Heu 


2.5.1  Parallel  polarization 


(2.9b) 


A  monochromatic  plane  wave  is  said  to  have  parallel  polarization  if  the  electric  held  is 
parallel  to  the  plane  of  incidence.  The  plane  of  incidence  is  dehned  by  the  wavevector  k 
and  the  normal  vector  to  the  surface  n;  i.e.,  k  and  n  lie  in  the  plane  of  incidence.  From 
Figure  2.1,  we  have  k  in  the  xy-plane  and  n  =  ±x,  thus,  the  plane  of  incidence  is  the 
xy-plane. 

Consider  a  parallel  polarized  incident  plane  wave  of  angular  frequency  u  propagating 
in  the  positive  x-direction.  Maxwell’s  equations  (2.1)  are  linear  PDFs,  thus,  the  total  wave 
inside  each  region  may  be  decomposed  into  rehected  and  transmitted  waves  with  the  following 
wavevectors: 

k^  =  +  ky^ij,  (2.10) 

where  kx,e  is  given  by  (2.8),  ky/  is  given  by  (2.7),  indicates  a  transmitted  wave  propagating 
in  the  +a;-direction,  and  “  indicates  a  rehected  wave  propagating  in  the  — x-direction;  notice 
that  there  is  no  rehected  wave  in  the  0th  region,  see  Figure  2.1.  The  magnetic  intensity  in 
each  region  is  given  by 

H^(r,t)  =  CiuE^exp  [i  (k^  ■  r  —  tat)]  z,  (2-11) 

where  E'^  is  the  complex  amplitude  associated  with  the  transmitted  wave,  Ef  is  the  complex 
amplitude  associated  with  the  rehected  wave,  and  Ef^^  =  0.  Substituting  (2.11)  into  (2.9a) 
yields 

E^(r,t)  =  E^exp  [i  (k^  ■  r  -  tat)]  [-ky^^xE  kx/y] .  (2.12) 
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From  (2.2b),  we  see  that  the  y-component  of  the  total  electric  held  and  the  total  magnetic 
intensity  are  continuous  across  the  interface.  It  is  convenient  to  dehne  a  new  symbol  for  the 
y-component  of  the  electric  held  evaluated  on  the  interface.  Let 

p 

xf  =  ±K,eEfexp  hs  ,  (2.13) 

s=i+l  . 

where  hi  is  the  thickness  of  the  £th  region  and,  for  convenience,  we  set  hi=Q  =  hi=p  =  0. 
In  (2.13),  xf=o,...,p-iJ  denotes  the  y-component  of  the  electric  held  at  the  interface  between 
regions  £  and  £  +  1  (the  interface  is  approached  from  the  £th  region),  and  x^=p  denotes  the 

^/-component  of  the  electric  held  at  the  interface  between  regions  p  and  p  —  1  (the  interface 

is  approached  from  region  p),  see  Figure  2.1.  Substituting  (2.11)  and  (2.12)  into  (2.2b),  and 
using  (2.13)  to  simplify  the  result,  yields 

^  =  Xt  +  Xi,  (2.14a) 

wi+i  =  We  {xj  -  Xi)  ,  (2.14b) 

for  £  =  0, . . . ,  p  —  1,  where 

we=^,  £  =  0,...,p.  (2.15) 

After  we  obtain  a  linear  system  for  the  perpendicular  polarization  case,  we  will  solve  the 
linear  system  given  by  (2.14),  see  Section  2.6. 

2.5.2  Perpendicular  polarization 

A  monochromatic  plane  wave  is  said  to  have  perpendicular  polarization  if  the  electric 
held  is  perpendicular  to  the  plane  of  incidence.  The  electric  held  in  each  region  is  given  by 

E^(r,  t)  =  exp  [i  ■  r  —  tat)]  z,  (2-16) 

where  is  given  by  (2.10),  and  the  ^  superscripts  have  the  same  meaning  as  in  Section  2.5.1. 
Also  as  in  Section  2.5.1,  we  set  =  0  because  there  is  no  rehected  wave  in  the  0th  region. 
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Substituting  (2.16)  into  (2.9b)  yields 


=  —  exp  [i  (k^  ■  r  -  cut)]  K/j]  ■  (2.17) 

Heu 

From  (2.2b),  we  see  that  both  the  total  electric  field  and  the  y-component  of  the  total 
magnetic  intensity  are  continuous  across  the  interface.  Let  the  electric  field  evaluated  on  the 
interface  be  denoted  by 


xf  =  Ef  exp 


Elk 


p 


■x,i  7  , 
s=e.+i 


h. 


(2.18) 


where  denotes  the  ;2-component  of  the  electric  field  at  the  interface  between  regions 

i  and  7+1  (the  interface  is  approached  from  the  7th  region)  and  xf=p  denotes  the  2:-component 
of  the  electric  field  at  the  interface  between  regions  p  and  p  —  1  (the  interface  is  approached 
from  region  p).  Substituting  (2.16)  and  (2.17)  into  (2.2b),  and  using  (2.18)  to  simplify  the 
result,  yields  (2.14),  where 


Wi  = - 7  =  0,...,p.  (2.19) 

Notice  that  the  linear  system  for  the  perpendicular  polarization  case  is  the  same  as  the  linear 
system  for  the  parallel  polarization  case,  but  the  definitions  of  xf  we  are  different. 


2.6  Linear  system 


The  traditional  approach  to  solving  the  linear  system  given  by  (2.14)  is  to  rewrite  it  as 

7  =  0,...,p-l,  (2.20a) 


'xi+i' 

=  Me 

'xT 

Xi+i_ 

Xi. 

where 


Me 


{we+i  +  We) 
(we+i  -  We)  ipe+i 


{we+i  -  We)  i)eli 
{we+i  +  We)  'ipe+i  ’ 


(2.20b) 


and  ipe  =  exp  {ikx,ehe)-  To  compute  xt i  we  iterate  (2.20a)  until  7  =  p  —  1  to  obtain 


xf/xt 

Xp  / Xo , 


Mp_iMp_2 


(2.21) 
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After  computing  xt  from  (2.21),  we  can  find  from  (2.20a).  The  approach  outlined  above 
is  the  standard  transfer  matrix  method,  but  unfortunately  it  is  numerically  unstable  because 
the  top  half  of  iVfr  grows  exponentially  and  the  bottom  half  of  iVfr  decreases  exponentially 
if  Im  [kx/h(\  7^  0.  To  avoid  the  numerical  instability,  we  must  reformulate  the  linear  system 
given  by  (2.14)  in  terms  of  ipi  or  alone.  If  Im  [kx,ehi]  is  large,  then  %jj£  may  cause  underflow 
errors  and  may  cause  overflow  errors.  Generally  speaking,  underflow  is  preferred  to 
overflow  because  when  underflow  occurs,  the  (normal)  number  is  rounded  to  the  nearest 
subnormal  number  or  to  0.0;  thus,  it  is  desirable  to  reformulate  the  linear  system  in  terms 
of  -ipi  instead  of  (see  Section  2.6.1). 

2.6.1  S'-matrix 

In  this  section,  we  present  a  particularly  simple  version  of  the  S'-matrix  formulation  of 
(2.14)  that  avoids  numerical  instabilities.  To  derive  the  S'-matrix,  we  write  a  scattering 
matrix  (S'-matrix)  for  an  “aggregate  layer”  consisting  of  0, ...  ,i  layers  to  obtain 


„(1,2)' 

^t 

'  0  ■ 

Xo. 

J2,l) 

ft 

J2,2) 

■^t 

xf. 

Using  (2.20)  to  eliminate  xf  from  (2.22)  and  comparing  the  result  to  (2.22)  with  i  ^  i  +  1 
yields 


J1.2)  _ 


we+1  -  We 


1  —  s 


(1,2) 


1  S 


(1,2) 


We+1  +  We 


1  —  s 


(1,2) 


1  S 


(1,2)' 


„(2,2)  _ 


2we+is); 


(2,2) 


We+1 


1  +  s 


(1,2) 


+  We 


Jl,2) 


-  -1 

z-Zli’e+ii 

(2.23a) 

=r  P^t+i  j 

(2.23b) 

for  £  =  0, ...  ,p  —  1,  where  =  0  and  =  1.  Substituting  i  =  p  into  (2.22)  yields 
Xo  /Xp  =  where  is  computed  recursively  from  (2.23b).  Using  (2.20)  to  compute  xf 


would  make  the  algorithm  numerically  unstable.  To  avoid  introducing  numerical  instability 
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in  the  computation  of  xf  i  we  eliminate  Xo  and  Xi+\  from  (2.20)  and  (2.22)  to  obtain 


Xj  = 


Wl^\ 


1  + 


+  We 


1  —  s 


(1,2) 


1  Xe+it 


(2.24a) 


for  £  =  p  —  1, . . . ,  0  and 

=  £  =  p,...,1.  (2.24b) 

Notice  that  Xe  only  depends  on  The  S'-matrix  algorithm  is  numerically  stable  because 

(2.23a)  and  (2.24)  only  depend  on 

Originally,  (2.23a)  and  (2.24)  were  derived  in  [16]  by  citing  the  general  scattering-theory 
paradigm  that  requires  existence  of  a  linear  relationship  between  x'e  and  i.e.,  x't  — 
I  and  then  substituting  it  directly  into  (2.20)  to  obtain  (2.23a)  and  (2.24a).  Arguably 
our  derivation  is  just  as  simple  as  in  [16]  but  follows  the  traditional  S'-matrix  formulation 
[15,  17]  more  closely. 

We  would  like  to  note  that  it  is  possible  to  formulate  an  S'-matrix  algorithm  where  Xe 
are  computed  directly  from  y+  [18,  19],  but  such  a  formulation  requires  recursive  com¬ 
putation  of  three  elements  of  an  S'-matrix  rather  than  just  one  element  in  our  formu¬ 
lation.  Moreover,  it  is  also  possible  to  obtain  formulas  that  directly  relate  xf  to  Xp 
from  our  formulation  by  simply  multiplying  out  (2.24),  i.e.,  xt  =  ^’^d 

Xe  =s}  ■■■4  ^)Xp,  where  / s\ 


2.7  Conserved  quantities 


In  the  case  of  the  RHM,  the  time-averaged  complex  Poynting  theorem  for  harmonic  fields 
is  given  by 

V  •  S  +  gfr)  +  +  2ia;  =  0,  (2.25a) 

where  S  =  x  H*  is  the  complex  Poynting  vector  and 
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«(")  =^E-E*  =  ^  ||E||% 

4  4  "  "  ’ 

(2.25b) 

(2.25c) 

Q''’=^E.E*  =  ^A||Ef, 

(2.25d) 

IlHf. 

(2.25e) 

In  (2.25),  is  the  real  time-averaged  electric  density,  is  the  real  time-averaged  mag¬ 
netic  density,  and  represent  time-averaged  electric  and  magnetic  losses,  respectively 
(e.g..  Joule  heating  [37,  Sec.  2.19,  Sec.  2.20]).  Substituting  the  total  electric  held  and  the 
total  magnetic  intensity  into  (2.25b)  and  (2.25c),  respectively,  yields 

uP  =1  (ijE+lt  +  ||E,-f  +  2Re  [E+  •  E,-*])  ,  (2.26a) 

'4”>  =d  (||H+|f  +  ||H,-  f  +  2Re  [H+  ■  H,-*])  ,  (2,26b) 

where  Re  denotes  the  real  part. 

In  the  case  of  the  LHM,  the  complex  Poynting  theorem  for  harmonic  helds  given  by  (2.25) 
is  mathematically  correct.  However,  the  identihcation  of  the  real  electric  density  (2.25b) 
and  the  real  magnetic  density  (2.25c)  is  troublesome  because  both  are  negative.  It  was 
pointed  out  by  Veselago  [20]  that  the  LHM  must  be  accompanied  by  frequency  dispersion, 
in  which  case  the  real  electric  density  and  the  real  magnetic  density  are  not  given  by  (2.25b) 
and  (2.25c),  respectively.  Moreover,  simultaneously  negative  permittivity  and  permeability 
occur  very  near  resonance  and  there  is  therefore  no  frequency  interval  for  the  LHM  where 
permittivity  and  permeability  may  be  reasonably  approximated  by  a  constant.  For  a  more 
detailed  discussion  see  [23,  38,  39]. 

Another  conserved  quantity  is  the  fundamental  invariant  in  multilayers  (FIM)  [40,  41], 
given  by 

Wi+i  {i^e+ixt+iT  -  =we  (xiT  -  (xf)^  ,  (2-27) 
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for  £  =  0, ...  ,p  —  1.  The  FIM  is  a  product  of  the  continuity  conditions  for  the  electric  held 
(2.14a)  and  magnetic  intensity  (2.14b).  However,  the  FIM  is  not  an  energy  conservation 
statement  because  it  contains  and  instead  of  \xt'^  and  In  our  view, 

the  FIM  is  particularly  interesting  because  its  structure  is  similar  to  that  of  the  space-time 
interval  of  special  relativity,  ds^  =  dx^  —  c^dt^,  where  c  is  the  speed  of  light.  Moreover,  it 
has  been  pointed  out  in  [42]  that  many  results  associated  with  wave  propagation  through 
planar  stratihed  media  are  more  easily  derived  through  an  analogy  with  special  relativity. 
In  this  paper,  we  don’t  pursue  the  analogy  between  wave  propagation  though  a  multilayer 
stack  and  the  theory  of  special  relativity  any  further,  but  we  do  want  to  stress  that  this 
analogy  is  not  a  mere  coincidence. 

2.7.1  Energy  densities  for  parallel  polarization 

It  is  convenient  to  introduce  a  new  symbol  for  the  transverse  component  (the  y-component) 
of  the  electric  field  as  a  function  of  distance,  x,  into  the  multilayer  stack.  For  £  =  0, . . .  ,p. 


F^(a;)  =  [±ik^^(,x\ 


(2.28) 


then. 

Re  [F+(a;)F7*(a;)]  =  -  |A;.y|'Re  . 

Substituting  (2.12)  into  (2.26a)  and  using  (2.29)  to  simplify  the  result  yields 


(2.29) 


r+(i)f +  |r,-(u 


+2  1 


Re  [r+(i)r, -•(!)]  .  (2.30) 


Substituting  (2.11)  into  (2.26b)  and  using  (2.29)  to  simplify  the  result  yields 


^  (|r+(i)|"  +  |r,-(i)p  -  2Re  [r,+(i)r7-(i)]) ,  (2.31) 
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where  wg  is  given  by  (2.15). 


2.7.2  Energy  densities  for  perpendicular  polarization 

Again,  it  is  convenient  to  introduce  a  new  symbol  for  the  transverse  component  (the 
;2-component)  of  the  electric  held  as  a  function  of  distance,  x,  into  the  multilayer  stack.  For 
£  =  0, . . .  ,p,  let 

r^(a;)  =  exp  [±ikx/x]  (2.32) 

then, 

|r^(a^)r  =  \Ef\‘^  exp  {T2lm[kx,i\x), 

Re  [r+(a;)r7*(a;)]  =Re  . 

Substituting  (2.16)  into  (2.26a)  and  using  (2.33)  to  simplify  the  result  yields 

4\x)  =  I  [|r+(i)|"  +  [r,-(i)|"  +  2Re  [r+(i)r,--(i)]] .  (2.34) 

Substituting  (2.17)  into  (2.26b)  and  using  (2.33)  to  simplify  the  result  yields 

=  +  (|ryx)t  +  |r,-wp) 

r  "  i^)  W]  ,  (2.35) 

where  wg  is  given  by  (2.19). 

2.8  Transmission  and  reflection  coefflcients 

The  transmission  coefficient,  T,  and  the  reflection  coefficient,  R,  are  given  by 

y  ^  Re[S+]  ■  X 
Re[S+]  ■  x’ 

Re[S;]-x 
Re[S+]  ■  x’ 

with 


(2.36a) 

(2.36b) 
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S+  =  iE+  X  Hr  and  Sj  =  X  Hf, 

where  it  is  understood  that  and  *  are  evaluated  at  the  interface  between  regions  p 
and  p—1  (the  interface  is  approached  from  region  p),  and  Eg  and  Hg  *  are  evaluated  at  the 
interface  between  regions  1  and  0  (the  interface  is  approached  from  the  0th  region). 

In  the  case  of  the  parallel  polarization  state,  substituting  (2.11)  and  (2.12)  into  (2.36), 
and  using  (2.13)  to  simplify  the  result,  yields 


T  = 


R  = 


I  ^x,Q 


i^x,0_ 

Xo 

|2 

o| 

Xp" 

JL 

2 

.+ 

p 

(2.37a) 

(2.37b) 


In  the  case  of  the  perpendicular  polarization  state,  substituting  (2.17)  and  (2.16)  into 
(2.36),  and  using  (2.18)  to  simplify  the  result,  yields 


T  =  ^Re 


R  = 


Xr 


xh 


Xi 


Xi 


(2.38a) 

(2.38b) 


The  transmission  and  reflection  coefficients,  given  by  (2.37)  for  the  parallel  polarization 
state  and  by  (2.38)  for  the  perpendicular  polarization  state,  are  valid  for  both  a  right-  and 
a  left-handed  material. 


2.9  Multilayer  classes 

Python  is  a  multi-paradigm  programing  language  that  supports  object-oriented  program¬ 
ing,  structured  programing,  and  a  subset  of  functional  and  aspect-oriented  programing  styles. 
There  is  a  large  number  of  numerical  libraries  available  for  use  with  Python.  We  chose  to 
use  a  numerical  library  called  SciPy  [43]  for  numerical  computations  because,  in  our  opinion, 
a  reader  familiar  with  MATLAB"*"^  and/or  Fortran  90/95  will  hnd  SciPy  a  very  natural  and 
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easy-to-use  library. 


In  order  for  our  multilayer  classes,  namely  Boundary  and  Layer,  which  are  collectively 
called  openTMM,^  to  be  as  useful  as  possible  to  the  scientihc  community,  we  paid  particular 
attention  to  the  readability,  usability,  and  maintainability  of  the  code.  Both  classes  are 
implemented  in  an  object-oriented  programing  style  as  described  below. 

The  Boundary  class  is  meant  to  be  a  base  class  {superclass  in  the  Python  lexicon)  that 
will  be  inherited  by  the  derived  classes  {subclasses  in  the  Python  lexicon).  The  derived 
classes  perform  “high-level”  computations  such  as  computing  the  energy  density  and  the 
transmission  and  reflection  coefficients.  The  derived  Layer  class  inherits  the  Boundary  and 
computes  the  quantities  described  in  Table  2.1.  The  beneht  of  using  inheritance  in  our 
multilayer  calculations  is  that  other  developers  may  extend  the  Layer  class  or  write  their 
own  derived  class  to  compute  the  desired  quantity  of  interest  without  having  to  implement  the 
low-level  code,  e.g.,  the  code  for  computing  kx/  and  the  S'-matrix.  The  Boundary  superclass 


Table  2.1;  The  first  column  contains  the  name  (as  it  appears  in  the  code)  of  the  object 
attribute  (method)  of  the  class  Layer,  the  second  column  contains  a  description  of  the 
method,  and  the  third  column  contains  references  to  the  section  where  a  more  detailed 
description  may  be  found. 


Name 

Description 

Refs. 

field 

Transverse  component  of  the  electric  held  as  a  function 
of  distance,  r^(a;) 

2.7.1 

2.7.2 

energy 

Electric/magnetic  energy  density  as  a  function  of 
distance, 

2.7.1 

2.7.2 

loss 

Electric/magnetic  losses  as  a  function  of  distance, 

2.7 

divPoynting 

Divergence  of  the  Poynting  vector  as  a  function  of 
distance,  V  ■  S(a;) 

2.7 

FIM 

EIM  at  each  boundary  interface 

2.7 

FIMvsDist 

FIM  as  a  function  of  distance 

2.7 

TRvsFreq 

TRvsAngle 

TRvsFreqAndAngle 

Transmission  and  rehection  coefficients  as  a  function  of 
frequency  /  =  u/2n  and/or  angle  of  incidence  0,  i.e., 
{T{f),R{f)},  {T(0),i?(0)},  {T(/,0),i?(/,0)} 

2.8 

^openTMM  is  an  open-source  software  distributed  under  the  MIT  license  and  is  available  from  http://pypi. 
python . org/pypi/ openTMM. 
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computes  a  “minimal”  set  of  “basic”  quantities,  see  Table  2.2,  that  are  used  by  the  Layer 
subclass.  Each  function/method  in  the  Boundary  and  Layer  class  contains  a  documentation 
header  {docstring  in  the  Python  lexicon),  which  describes  the  function/method  in  detail 
and  includes  an  example  of  its  use.  To  access  the  docstrings,  the  user  may  use  Python’s 
help  function  or  if  more  user  friendly  formatting  is  desired,  the  user  may  use  SciPy’s  info 
function.  For  example,  the  docsting  for  Layer .  energy  function  may  be  accessed  via 

»>  help (openTMM. Layer. energy) 

»>  scipy .  info  (openTMM .  Layer .  energy) 

and  all  docstings  contained  in  a  class  may  be  accessed  via 

»>  help  (openTMM.  ClassName) 

»>  scipy .  inf o (openTMM.  ClassName) 

where  ClassName  is  either  Boundary  or  Layer.  This  interactive  documentation  feature  of 
Python  makes  it  a  very  convenient  language  to  use  and  largely  eliminates  the  need  to  produce 
separate  code  documentation.  The  help/scipy .  inf  o  functions  are  similar  to  the  Manual 


Table  2.2:  The  first  column  contains  the  name  (as  it  appears  in  the  code)  of  the  object 
attribute  of  the  class  Boundary,  the  second  column  contains  a  description  of  the  attribute, 
and  the  third  column  contains  references  to  a  section  and/or  equation  where  a  more  detailed 
description  of  the  attribute  may  be  found. 


Name 

Description 

Refs. 

self  .h 

Thickness  of  each  layer, 

(2.13),  (2.18) 

self . epsRel 

Relative  permittivity  of  each  region,  e^/cvacuum 

Section  2.4 

self .muRel 

Relative  permeability  of  each  region,  /i^Z/Xvacuum 

Section  2.4 

self .pol 

Polarization  state 

Section  2.5 

self .kx 

x-component  of  the  wavevector, 

(2.8) 

self .  w 

Scaled  self  .kx  (polarization  dependent), 

(2.15),  (2.19) 

self . chiPlus 

Transverse  component  of  the  electric  held  evaluated 
on  the  interface,  xt /Xp 

(2.13),  (2.18) 

self . chiMinus 

Transverse  component  of  the  electric  held  evaluated 
on  the  interface,  Xe  /Xp 

(2.13),  (2.18) 
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pager  utils  (man  pages)  of  Unix-like  operating  systems;  could  one  imagine  using  a  Linux 
shell  without  man  python? 

2.10  Python  and  nnmerical  efficiency 

There  is  some  concern  about  the  speed  of  computations  in  Python  because  it  is  byte- 
compiled,  not  a  compiled  language  such  as  Fortran  90/95  or  C/C-I--I-.  However,  in  our 
opinion,  the  code  readability  (less  error-prone  syntax),  flexibility  (effortless  integration  with 
other  software)  and  ease-of-use  of  Python  (leading  to  shorter  development  times)  in  many 
cases  outweigh  any  performance  benefits  of  compiled  languages.  An  interested  reader  may 
consult  [44-47]  for  a  fuller  discussion  of  why  Python  is  a  language  of  choice  for  scientific  soft¬ 
ware  development.  Typically,  computationally  intensive  routines  in  Python  are  implemented 
in  compiled  languages  and  therefore,  the  difference  in  computation  time  between  Python  and 
complied  languages  is  acceptable  for  many  applications  [45-48].  In  the  Python  lexicon,  the 
mixing  of  programing  languages  is  called  the  Pythonic  approach]  this  is  the  approach  we  use 
with  the  computationally  intensive  part  of  the  Boundary  superclass. 

It  is  relatively  obvious  that  the  computationally  intensive  part  of  the  Boundary  superclass 
is  the  computation  of  xf:  i-®-)  the  solution  of  the  linear  system  described  in  Section  2.6. 
Therefore,  the  computation  of  xf  is  implemented  in  Fortran  90  and  the  Python  bindings 
are  built  by  F2PY  [49]  (F2PY  is  now  part  of  SciPy).  However,  implementing  “workhorse 
functions”  in  a  compiled  language  reduces  the  readability  and  maintainability  of  code  to 
some  extent.  Therefore,  we  strongly  encourage  developers  to  only  implement  workhorse 
functions  in  compiled  languages  when  they  lead  to  severe  bottlenecks.  It  is  often  the  case  that 
bottlenecks  can  only  be  identified  after  code  profiling  (performance  analysis).  For  example, 
it  is  not  obvious  that  the  square  root  function  in  the  computation  of  is  relatively  time- 
consuming.  The  computation  of  kx,e  is  relatively  expensive  because  SciPy’s  square  root 
function,  scipy.sqrt,  does  an  element-by-element  analysis  of  the  input  array  to  find  if  it 
contains  any  real  elements  less  than  zero.  If  a  real,  less-than-zero  element  is  found,  SciPy 
converts  the  whole  input  array  to  a  complex  data  type  and  passes  it  to  NumPy  [50],  which 
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uses  an  efficient  C  code  to  compute  the  square  root.  In  our  case,  SciPy’s  time-consuming 
element-by-element  analysis  is  unnecessary  because  of  a  priori  knowledge  about  kx^e.-,  see 
Section  2.4.  We  could  avoid  scipy.sqrt  by  directly  using  NumPy’s  square  root  function, 
but  this  is  not  the  most  convenient  approach  because  NumPy’s  square  root  function  of 
a  complex  number  =  \z\e^^  returns  where  — tt  <  6  <  tt,  but  (2.8)  requires 

that  \m[kx^e\  >  0.  To  avoid  this  inconvenience,  we  choose  to  implement  our  own  square 
root  subroutine,  cmplx_sqrt,  which  returns  the  square  root  in  an  appropriate  quadrant  as 
required  by  (2.8).  The  cmplx_sqrt  is  implemented  in  Fortran  90  with  Python  binding  build 
by  F2PY  and  depends  on  Fortran’s  intrinsic  square  root  function,  SQRT. 

To  conhrm  that  the  run-time  of  the  Python  Boundary  superclass  is  acceptable,  we  com¬ 
pared  it  to  a  Boundary  class  implemented  in  pure  Fortran  90.  From  Figure  2.2,  we  see  that 
for  a  large  number  of  layers  (>  300)  the  Python  code  is  only  25  percent  slower  than  the 
pure  Fortran  90  code.  However,  for  a  small  number  of  layers  (<  20)  the  Python  code  is 
about  10  times  slower  than  the  pure  Fortran  90  code,  see  inset  in  Figure  2.2,  but  this  is 
not  major  concern  because  such  a  small  number  of  layers  has  an  absolute  execution  time 
about  a  second  or  so  in  Python.  We  believe  that  the  run-time  discrepancy  between  a  small 
and  a  large  number  of  layers  is  caused  by  SciPy’s  overhead  cost,  which  does  not  increase 
signihcantly  as  a  function  of  array  size. 

2.11  Numerical  stability  and  accuracy 

To  demonstrate  the  numerical  stability  and  accuracy  of  openTMM,  we  numerically  checked 
the  complex  Poynting  theorem  given  by  (2.25a),  the  fundamental  invariant  in  multilayers 
given  by  (2.27),  and  the  de  Hoop  reciprocity  theorem  [51,  Section  6].  The  de  Hoop  reciprocity 
theorem  states  that  if  eo  =  and  Hq  =  /ip,  then  the  transmitted  wave  is  unaffected  by  a 
180  degree  rotation  of  the  multilayer  stack  around  the  2:-axis,  see  Figure  2.1.  We  measure 
the  accuracy  of  a  computed  quantity  in  terms  of  the  number  of  significant  digits  it  agrees 
with  the  theoretical  value  and  we  denote  this  measure  of  accuracy  by  (5„.  Approximately,  5^ 


number  of  layers 

Figure  2.2:  The  ratio  of  total  computational  time  required  to  compute  T(/j,  0^)  and  R{fi,  0j), 
where  1  <  {i,j}  <  10^,  using  openTMM  and  a  pure  Fortran  90  code.  Each  multilayer  stack  is 
composed  of  the  same  number  of  pseudorandom  layers  of  the  following  types:  right-handed 
layers  with/without  absorption  and  left-handed  layers  with/without  absorption. 

is  given  by 

X  .  r  ,  -  ^Re  ,  v-vi^  \ 

dy  ^  mm  <  —  log  -  ,  —  log  -  >  (2.39) 

1^  u  J 

where  v  is  the  theoretical  value  and  {uRe,him}  are  the  numerically  computed  value.  For  a 
numerical  check  of  the  complex  Poynting  theorem  URe  is  given  as  —Re  [V  ■  S]  / 
and  him  is  given  as  — Im  [V  ■  S]  /  [2a;  For  the  FIM  test,  URe  (uim)  is  the  ratio 

of  the  the  real  (imaginary)  part  of  the  left-hand  side  to  the  real  (imaginary)  part  of  the 
right-hand  side  of  (2.27).  For  the  de  Hoop  reciprocity  test,  URe  =  Re  [xg]  /Re  \Xp'\  and 
him  =  Im  [xo]  /Im  \Xp'\i  where  Xo  i^  fhe  transmitted  wave  before  the  180  degree  rotation 
of  the  multilayer  stack  and  Xp  is  the  transmitted  wave  after  the  rotation.  For  all  three 
numerical  checks,  v  =  1  and  all  computations  are  performed  in  double-precision  (^  16 
signihcant  digits).  From  Figure  2.3,  we  see  that  the  three  numerical  checks  are  satished  with 
an  accuracy  of  >  12.  Despite  the  fact  that  some  of  the  layers  in  the  stack  chosen  for 
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Figure  2.3:  6^  is  shown  for  a  multilayer  stack  composed  of  the  same  number  of  pseudorandom 
layers  of  the  following  types:  right-handed  layers  with/without  absorption  and  left-handed 
layers  with/without  absorption. 

Figure  2.3  have  very  high  absorption,  Im  [kx,ehi]  ~  30,  we  see  that  6^  does  not  decrease  as 
a  function  of  distance  into  the  stack,  which  conhrms  that  our  S'-matrix  algorithm  is  indeed 
numerically  stable.  The  composition  of  the  multilayer  stack  is  summarized  in  Table  2.3.  To 
produce  Figure  2.3,  we  used  a  normally  incident  plane  wave  of  frequency  100  GHz  and  the 
hrst  pathological  case  scheme,  see  Section  2.4.2. 

2.12  Conclusions 

A  numerically  stable  S'-matrix  algorithm  for  electromagnetic  wave  propagation  through 
planar  stratihed  media  composed  of  a  right-handed  and/or  left-handed  material  has  been 
implemented  in  Python.  Pathological  cases  caused  by  an  unphysical  approximation  of  zero 
absorption  have  been  carefully  examined  and  numerically  circumvented  (see  Section  2.4.2). 
The  numerical  computations  were  implemented  in  an  object-oriented  programing  style  by 
dividing  them  into  two  classes.  Boundary  and  Layer.  The  Boundary  class  performs  com- 
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Table  2.3:  The  height,  h,  relative  permittivity,  erei,  and  the  relative  permeability,  /irei,  of 
each  layer  were  pseudorandomly  chosen  from  the  intervals  shown  in  the  table.  E.G.,  from 
the  second  line  of  the  table,  we  see  that  125  layers  have  thickness  between  1  mm  and  10  mm, 
and  relative  permittivity/permeability  between  —10  and  —1.  A  500-layer  stack  was  used 
for  the  Poynting  theorem  and  the  FIM  test.  For  the  de  Hoop  reciprocity  theorem  test,  a 
multilayer  stack  consisting  of  5,  9,13,...,  501  layers  was  used. 


Test 

of  layers 

h  (mm) 

<61 

f” 

^rel 

h-rel 

/^rel 

Poynting  Thm 
&  FIM 

125 

[1,10] 

[1, 10] 

0 

11.10] 

0 

125 

[1, 10] 

[-10,-1] 

0 

[-10,-1] 

0 

125 

[1,10] 

[-10,-1] 

[0.01,0.1] 

[-10,-1] 

0 

115 

[1,10] 

[1, 10] 

0 

[1,10] 

[0.01,0.1] 

10 

[10, 15] 

[2, 10] 

|0.2| 

[1,10] 

0 

de  Hoop  Thm 

{!,...,  125} 

[1,10] 

[1, 10] 

0 

[1,10] 

0 

{!,...,  125} 

[1,10] 

[-10,-1] 

0 

[-10,-1] 

0 

{!,...,  125} 

[1,10] 

[-10,-1] 

[0.01,0.1] 

[-10,-1] 

0 

{!,...,  125} 

[1,10] 

[1, 10] 

0 

[1,10] 

[0.01,0.1] 

1 

[1,10] 

[2, 10] 

[0.1] 

[1,10] 

0 

putationally  intensive  calculations,  namely  the  solution  of  the  linear  system  described  in 
Section  2.6.1  and  the  square  root  of  The  workhorse  functions  of  the  Boundary  class 
were  implemented  in  Fortran  90  in  order  to  avoid  computational  bottlenecks.  The  Layer 
class  performs  high-level  calculations,  such  as  calculation  of  r=*"(a;),  and 

FIM.  The  code  has  been  tested  and  is  accurate  to  ~  12  signihcant  digits  (see  Section  2.11). 

We  hope  that  our  open-source  and  object-oriented  implementation  of  the  A-matrix  al¬ 
gorithm,  which  is  suitable  for  modern  applications  such  as  Anderson  localization  of  light 
and  perfect  lensing,  will  be  adopted  by  a  wide  scientihc  community.  At  the  very  least,  we 
hope  that  our  publicly  available  implementation  of  the  A-matrix  algorithm  will  encourage 
the  scientific  community  to  use  open-source  software,  thus  increasing  the  reproducibility  of 
scientihc  work. 
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3.1  Abstract 

Linear  response  laws  and  causality  (the  effect  cannot  precede  the  cause)  are  of  fundamen¬ 
tal  importance  in  physics.  In  the  context  of  classical  electrodynamics,  students  often  have  a 
difficult  time  grasping  these  concepts  because  the  physics  is  obscured  by  the  intermingling 
of  the  time  and  frequency  domains.  In  this  paper,  we  analyse  the  linear  response  laws  and 
causality  in  the  time  and  frequency  domains  with  the  aim  of  pedagogical  clarity.  We  will 
show  that  it  is  easy  to  violate  causality  in  the  frequency  domain  by  making  a  vanishing 
absorption  approximation.  Further,  we  will  show  that  there  can  be  subtle  differences  be¬ 
tween  Fourier  transforming  Maxwell  equations  and  using  a  monochromatic  source  function. 
We  discuss  how  these  concepts  can  be  obscured  and  offer  some  suggestions  to  improve  the 
situation. 

3.2  Introduction 

When  encountering  Maxwell’s  equations  in  matter  for  the  hrst  time,  students  are  faced 
with  many  conceptual,  as  well  as  mathematical,  difficulties.  In  the  time-domain,  the  four 
macroscopic  Maxwell’s  equations,  namely. 


V  ■  D  =  dvrp. 


V-B  =  0, 


(3.1) 
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are  straightforward  enough.  The  problem  arises  in  the  constitutive  relations  connecting  E 
and  D  or  B  and  H.  Whereas  the  time-domain  Maxwell  equations  are  real  and  involve 
only  real  quantities,  the  associated  response  functions  are  temporally  non-local  and  their 
very  dehnition  involves  integration  with  respect  to  time.  To  professionals,  this  issue  is  well 
understood  and  is  almost  always  side-stepped  notationally  by  writing  the  response  equations 
as 

D  =  eE  and  H  =  (3.2) 

and  mentally  juggling  the  time  and  frequency  domains.  Here,  for  simplicity,  we  have  assumed 
that  the  medium  is  linear,  isotropic  and  homogeneous  (LIH).  Experts  understand  that  (3.2) 
means  either 

(a)  Fourier-domain  relations  of  general  validity  (within  LIH  assumptions)  or 

(b)  time-domain  relations  valid  only  for  monochromatic  helds. 

But  to  simultaneously  present  (3.1)  and  (3.2)  to  students  (as  is  done  in  popular  textbooks 
[1-3])  can  mislead  and  confuse  them;  it  obscures  the  important  temporal  non-locality  of 
the  response  functions,  mixes  time  and  frequency  domain  concepts,  and  inserts  complex 
quantities  into  manifestly  real  equations.  To  avoid  a  possible  source  of  confusion,  we  use  the 
phrase  ‘Fourier-domain’  instead  of  ‘frequency-domain’  because  the  latter  may  also  refer  to  a 
time-domain  relation  with  a  monochromatic  source. 

Moreover,  students  fail  to  understand  the  relationship  between  absorption,  dispersion  and 
causality  (the  effect  cannot  precede  the  cause  [4]).  To  see  the  interplay  between  absorption, 
dispersion  and  causality  pictorially,  we  follow  Toll’s  ingenious  presentation  [5]  by  considering 
an  input  signal  S(t)  that  is  zero  for  t  <  0.  The  input  S(t)  is  a  weighted  sum  of  many 
sinusoidal  components,  such  as  Ci(t)  =  sin(a;it  +  ■^i),  each  of  which  extends  from  t  =  —oo 
to  t  =  oo.  Notice  that  the  weighted  sum  of  all  sinusoidal  components  produces  a  zero  input 
signal  for  t  <  0.  If  a  system  only  absorbs  one  component,  e.g.,  Ci(t),  without  affecting 
other  components,  then  the  output  of  such  a  system  would  simply  be  S(t)  —  Ci(t),  which  is 
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non-zero  for  t  <  0,  see  Figure  3.1.  Such  a  system  is  impossible  because  it  violates  causality; 
the  output  is  non-zero  before  the  onset  of  the  input  signal.  Therefore,  for  causal  systems, 
an  absorption  of  one  frequency  must  be  accompanied  by  phase  shifts  of  other  frequencies 
in  order  to  produce  a  zero  output  for  t  <  0,  and  the  necessary  phase  shifts  are  prescribed 
by  the  dispersion  relation.  Moreover,  the  converse  is  true  as  well;  namely,  a  phase  shift  of 
one  frequency  is  necessarily  accompanied  by  an  absorption  at  other  frequencies.  From  Toll’s 
argument,  we  conclude  that  it  is  impossible  to  design  a  physical  system  that  is  causal  and 
dispersionless.  Therefore,  when  one  speaks  of  dispersionless  media,  one  violates  a  sacred 
physical  principle  (causality). 


Figure  3.1;  The  input  signal,  S(t),  is  shown  along  with  the  output  signal,  S(t)  —  Ci(t),  of  a 
system  that  only  absorbs  the  Ci(t)  component  of  S(t),  without  affecting  other  components. 

This  paper  is  intended  not  only  for  instructors  who  teach  advanced  undergraduate  and 
beginning  graduate  students,  but  also  for  the  graduate  students  themselves.  The  mathe¬ 
matical  sophistication  needed  to  understand  this  paper  is  essentially  that  of  a  beginning 
graduate  student.  However,  throughout  the  paper  we  extensively  use  the  theory  of  the  tem¬ 
pered  distributions;  with  which  a  typical  beginning  graduate  student  may  not  be  familiar. 
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To  remedy  this  deficiency,  we  have  included  an  example-driven  tutorial  on  the  formal  theory 
of  tempered  distributions  in  the  appendix. 

3.3  Background 

The  time-domain  relationship  between  D  and  E  for  a  LIH  and  time-translationally  in¬ 
variant  medium  is  given  by® 

/OO 

x{t  —  t')E(t')  dt' ,  (3.3a) 

■OO 

where 


x{t  —  t')  =  a{t  —  t')Q{t  —  t')  and  — 


1,  t-t'>d 
0,  t-t'  <d 


(3.3b) 


The  appearance  of  the  Heaviside  step  function,  0,  in  the  definition  of  the  electric  suscep¬ 
tibility,  X,  reminds  us  that  the  polarization  vector,  P,  can  only  depend  on  the  past  values 
of  the  applied  electric  field,  E.  Therefore,  (3.3)  gives  a  causal  relationship  between  the 
displacement  field,  D,  and  the  applied  electric  field.  Taking  the  Fourier  transform  of  (3.3) 
yields 


D(a;)  =  E(a;)  -|-47rP(a;)  and  P(i^)  =  x('^)E('^))  (3-4) 


with  the  Fourier  transform  pair  given  by 


/(a;)  =  T[/(t)] 


f{t)  =  T-' 


1 


(3.5a) 

(3.5b) 


Under  suitable  mathematical  conditions,  the  principle  of  causality  in  the  time-domain 
translates  into  the  Kramers-Kronig  (KK)  relations  (the  Hilbert  transform  pair)  in  the 
Fourier-domain,  namely. 


®For  our  purposes,  it  is  more  convenient  to  work  with  (3.3)  than  with  D  =  e{t  —  dt' ,  where 

e{t  —  t')  =  d{t  —  t')  +  diTxit  —  t')- 
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Imx(a;)  =  ^  [Rex(r/)]  =  --/  dr/,  (3.6a) 

Rex(a;)  =  JC"'  [Imx(r/)]  =  dr/,  (3.6b) 

TtJ-oo 

where  Rex(a;)  and  Imx(a;)  are  the  real  and  imaginary  parts  of  x(ca),  respectively,  and 
denotes  the  Cauchy  principal  value  integral.  In  free  space,  by  definition,  x(t  —  f)  =  0  and 
consequently,  the  KK  relations  are  trivially  satisfied.  The  KK  relations  can  be  derived  in 
‘two  lines’  by  using  the  Fourier  representation  of  the  Heaviside  step  function  and  utilizing 
the  freedom  to  dehne  a{t  —  t')  for  t  —  t'<0[6].  A  lengthier  but  more  traditional  derivation  of 
the  KK  relations  is  available  in  [7].  This  derivation  is  self-contained  and  does  not  assume  a 
priori  knowledge  of  the  theory  of  functions  of  a  complex  variable.  For  a  historically  accurate 
account  on  how  the  KK  relations  were  hrst  derived,  see  [8]. 

In  general,  it  is  difficult  to  establish  an  equivalence  between  the  KK  relations  (3.6)  and 
causality  (3.3b).  If  x{^)  R  a  square-integrable  function,  then  the  Titchmarsch  theorem 
[9]  guarantees  that  (3.6)  and  (3.3b)  are  equivalent.  In  other  words,  x(t  —  t')  oc  0(f  —  t') 
if  and  only  if  x(a;)  satishes  (3.6).  The  square-integrability  requirement  on  x(a;)  may  be 
somewhat  relaxed.  It  can  be  shown  that  if  x{^)  R  bounded  but  E(t)  and  P(t)  are  square- 
integrable  functions,  then  (3.6)  and  (3.3b)  are  equivalent  [5].  Unfortunately,  these  conditions 
are  generally  not  satished,  in  fact,  they  are  not  even  satished  in  the  idealized  examples 
considered  in  this  paper.  Therefore,  we  must  part  with  our  naive  notion  that  y  is  a  classical 
function  and  treat  it  as  a  generalized  function  (distribution).  For  our  purposes,  it  will  be 
sufficient  to  treat  y(t)  as  a  tempered  distribution.  A  reader  not  familiar  with  the  formal 
theory  of  the  tempered  distributions  should  read  the  appendix  which,  for  our  purposes,  serves 
as  a  self-contained  tutorial  on  tempered  distributions.  For  the  reader’s  convenience.  Table  3.1 
provides  a  summary  of  the  notation  introduced  in  the  appendix.  From  this  point  on,  unless 
explicitly  noted  otherwise,  the  symbol  y(t)  should  be  interpreted  as  a  tempered  distribution, 
i.e.  x{t)  ^  S'.  Consequently,  the  convolution  integral  in  (3.3a)  and  the  Fourier  transform 
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Table  3.1:  A  brief  description  of  the  symbols  introduced  in  the  appendix. 


Symbol 

Description 

The  Schwartz  class 

S' 

The  class  of  tempered  distributions 

(/,0)  =  /!L/(f)0(f)  dt 

Tempered  distribution  f]fES'  and  0  G  S' 

(/  *  ^)  {t)  =  JZo  fit  -  t')'^it')  dF 

Convolution  of  /  with  f  ^  S'  and  -ip  E  S 

supp  [/] 

Support  of  f]  f  E  S' 

pair  given  by  (3.5)  should  also  be  interpreted  in  a  distributional  sense;  see  appendix  A. 1.1. 
Strictly  speaking,  the  Hilbert  transform  pair  given  by  (3.6)  should  also  be  interpreted  in  a 
distribution  sense.  However,  for  our  purposes,  it  will  be  sufficient  to  think  of  the  Hilbert 
transform  as  a  ‘generalized’  convolution  of  the  singular  function  l/ta  (Hilbert  kernel)  with 
a  tempered  distribution  f{u).  For  a  more  mathematical  treatment  of  the  Hilbert  transform 
of  the  tempered  distributions,  see  [10-13].  In  1958,  Taylor  [14]  (also,  see  discussion  in  [15]) 
rigorously  established  the  equivalence  between  causality  and  the  KK  relations  when  x{^)  is 
a  tempered  distribution.  This  is  a  very  important  result  that  we  will  use  repeatedly.  We  will 
force  x('^)  ^  S'  to  satisfy  the  KK  relations,  thereby  guaranteeing  that  x(t  —  t')  vanishes  for 
t  —  t'  <  0.  In  other  words,  Taylor’s  result  gives  causality  meaning  in  the  Fourier-domain. 

To  explore  these  fundamental  issues  in  a  pedagogical  context,  we  will  consider  two  text¬ 
book  models  that  are  used  to  derive  or  explain  the  dispersion  of  electromagnetic  waves. 
The  first  is  the  Drude  model  of  conduction  in  metals.  The  second  is  the  damped  harmonic 
oscillator,  which  arises  in  semi-classical  models  of  atomic  absorption.  The  former  is  just  a 
damped  harmonic  oscillator  model  with  a  spring  constant  of  zero,  and  where  we  interpret 
the  damping  in  terms  of  electron  collisions. 
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3.4  Drude  model 


In  a  metal,  the  motion  of  a  conduction  electron  of  charge  — e  and  mass  m  under  the 
influence  of  an  electric  held  E  is  given  by 

d^r  1  dr 


- = - E 

df2  T  dt  m  ’ 


(3.7a) 


where  r  is  the  collision  mean  free  time  [16].  Substituting  S{t  —  t')n  for  E  and  taking  the 
Fourier  transform,  see  appendix  A.  1.1,  yields 


u  {  u  -\ —  )  giu)  =  — E(a;)  and  E(a;)  =  e 
T  m 


n, 


(3.7b) 


where  g(cj)  denotes  the  Green’s  function,  n  =  ^  (x  +  y  +  z),  and  {x,  y,  z}  is  the  standard 
basis  for  the  three-dimensional  space.  Using  P(a;)  =  —nceg{u),  we  see  that  the  electric 
susceptibility  is  given  by 


u(u+-]x{co)  =  -^^,  (3.8a) 

\  T  J  m 

where  ric  denotes  the  conduction  electron  density  (assumed  to  be  constant).  When  solving 
(3.8a)  for  x(a;)  we  must  remember  that  x{^)  is  a  tempered  distribution;  in  other  words,  we 
seek  a  distributional  solution  to  (3.8a).  The  distributional  solution  consists  of  two  parts;  the 
particular  part,  Xpi^)^  given  by  the  ‘naive  division’  of  (3.8a),  namely. 


where  do 
satishes 


Xp(^^)  = 


O-Q 


(3.8b) 


u[tu  + 1) 

e’^UcT/m  is  the  static  conductivity,  and  the  homogeneous  part,  Xh(i^),  which 


Xh(w),  0(c^)^  =  (0, 0(cn))  (3.8c) 

for  all  (f)  E  S.  The  solution  to  (3.8c)  is  of  the  form  of  (A. 28)  and  not  of  (A. 36)  because 
u  d-i/r  never  equals  zero  for  real  u  and  hnite  r.  In  other  words,  supp  [Xh(i^)]  =  {0}  and  not 
{0,  — i/r}.  Therefore,  as  shown  in  the  hrst  example  in  appendix  A. 1.2,  the  form  of  the  full 
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solution  to  (3.8a)  is  given  by 


X{u)  =  -no  — - —  +  Xh{uj)  ,  where  Xh{uj)  =  coS{u).  (3.8d) 

\a;(ra;  +  i)  / 

To  find  the  unknown  constant  Cq,  we  require  that  the  response  of  the  system  be  causal. 
According  to  Taylor  [14],  the  causality  requirement  is  the  same  thing  as  requiring  that 
Imy  =  TC  [Re)^,  which  yields  cq  =  —n.  Therefore,  the  causal  electric  susceptibility  is  given 
by 

X{u)  =  -no  (  ^  -7r(5(a;)  )  ,  (3.9a) 

X{t  -  t')  =  [x(a;)]  =  uo  (^1  -  e“^  j  0(t  -  t').  (3.9b) 

Notice  that  x(t  —  t')  is  zero  for  t  <t'  and  increases  monotonically  to  its  maximum  value  of 
(To  for  t  >  t' .  The  rate  of  the  increase  is  controlled  by  the  collision  mean  free  time.  In  the 
next  two  subsections,  we  will  study  the  effects  of  x(t  —  t')  on  the  response  functions  when  a 
constant  or  monochromatic  electric  field  is  applied. 

3.4.1  A  constant  electric  field 

Suppose  now  that  a  constant  electric  field,  Eq,  was  turned  on  at  =  0.  Substituting 
E(T)  =  Eo0(T)  and  (3.9b)  into  (3.3a)  yields 

P{t)  =  (To  —  r  +  re“^  j  Eo0(t),  (3.10a) 

and  the  current  density  is  given  by 

J(t)  =  —  =  (To  (^1  -  e“- j  Eo0(t).  (3.10b) 

The  time  dependence  in  (3.10b)  reminds  us  that  when  we  connect  a  wire  to  a  battery,  the 
current  in  the  wire  does  not  instantaneously  reach  its  Ohm’s  law  steady-state  value.  The 
rate  at  which  the  current  reaches  the  steady-state  value  is  controlled  by  the  collision  mean 
free  time  r,  as  evident  from  the  exponential  term  in  (3.10b).  Of  course,  (3.10b)  reduces  to 
Ohm’s  law,  J  =  (ToEq,  when  t  3>  r. 
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3.4.2  A  monochromatic  field 


Another  simple  situation  is  a  monochromatic  driving  field  of  angular  frequency  cUd  that 
has  existed  since  the  beginning  of  time  (t'  =  —oo).  Substituting  E(f')  =  Eocos(a;dt')  and 


(3.9b)  into  (3.3a)  yields 


D(i)  =  AiEq  cos  (cudt)  +  A2E0  sin  {udt) , 


(3.11a) 


where 


Ai  —  1 


dTrUoT 
l  +  r2a;2 


and  Ao  = 


(3.11b) 


1,22  ^  /I  I  2  2^  • 

1  +  cud  (1  + 

From  (3.11a),  we  see  that  D(t)  has  one  component  that  oscillates  in-phase  and  one  com¬ 
ponent  that  oscillates  out-of-phase  with  the  applied  field.  The  out-of-phase  oscillations  are 
caused  by  collisions  of  the  electrons  with  the  ions  (absorption).  We  can  put  (3.11)  into  a 
form  of  (3.2)  if  we  let 


e  =  Ai  -I-  1^2  and  E  =  Eoe 


(3.12a) 


then,  (3.11)  becomes 


D(t)  =  Re  ( e  E 


(3.12b) 


Notice  that,  in  general,  D(t)  7^  Re(e)Re(E).  Thus,  what  professionals  mean  by  (3.2)  is  really 
(3.12)  if  they  are  talking  about  a  monochromatic  applied  field.  But  what  is  e?  To  answer 
this  question,  let’s  rewrite  (3.12a)  in  an  illuminating  form,  namely. 


e  —  1  _  ^  _  (Jo 

Ud{TUd  +  iy 


(3.13) 


Notice  that  (3.13)  has  a  simple  pole  at  cud  =  0.  The  root  cause  of  this  pole  is  our  assumption 
that  the  medium  is  homogeneous,  and  the  definition  of  homogeneity  clearly  depends  on 
the  spatial  wavelength  of  the  interrogating  field  [17,  18].  Thus,  (3.13)  may  be  devoid  of 
physical  reality  near  cud  =  0.  Be  that  as  it  may,  the  resemblance  between  xi^d)  and  x(a;)  is 
uncanny  if  we  replace  with  u  in  (3.13).  But  is  such  a  comparison  of  x  and  x  meaningful? 
Strictly  speaking,  it  is  not,  because  y  is  a  classical  (ordinary)  function  and  y  is  a  tempered 
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distribution  (generalized  function).  Nevertheless,  we  want  to  compare  these  two  objects  by 
‘promoting’  y  to  be  a  generalized  function.  Let  X{u)  =  ~  then 


{X(cn),0(a;)) 


7iaoS{u)(l){u)  du  =  7r(To0(O) 


(3.14) 


for  all  4)  E  S.  From  (3.14),  we  see  that  supp[X(a;)]  =  {0},  i.e.  the  ‘promoted’  x{^)  and 
x(a;)  differ  in  a  neighbourhood  of  a;  =  0.  In  the  laboratory,  we  actually  make  a  Fourier- 
domain  ‘measurement’  in  the  time-domain  by  driving  the  system  for  a  very  long  time  with  a 
monochromatic  held  [19].  Thus,  we  experimentally  measure  ^(cad)  and  not  x{^)-  Moreover, 
it  is  the  x{uj)  that  satishes  the  KK  relations  and  not  the  x(co’d),  even  if  we  ‘promote’  x(a;d) 
to  a  generalized  function.  This  crucial  difference  between  x{u)  and  x(a;d)  is  often  missed  by 
students. 

A  mathematically  inclined  reader  might  object  to  the  usage  of  E(t')  =  Eq  cos(a;d^0  for  a 
driving  held.  Clearly,  cos{udt')  ^  S  but  we  have  only  dehned  the  convolution  (A.  18)  between 
a  distribution  in  S'  and  a  function  in  S.  Of  course,  physically  we  know  that  a  starving  grad¬ 
uate  student  had  to  turn  the  held  on  and  oh.  If  we  assume  that  the  held  was  turned  on/oh  in 
a  smooth  and  slow  fashion,  then  it  could  be  approximated  by  cos(a;dC)  exp  [—  (at')^] ,  which 
does  belong  to  the  space  S.  A  similar  physical  argument  can  be  used  to  justify  the  usage  of 
the  Heaviside  step  function  in  Section  3.4.1.  Moreover,  under  certain  conditions,  it  is  pos¬ 
sible  to  dehne  convolution  of  two  tempered  distributions  [20].  This  dehnition  is  beyond  the 
scope  of  this  paper  but,  if  used,  it  would  eliminate  the  need  for  the  above  physical  argument. 


3.5  Plasma 


A  limiting  case  of  the  Drude  model  is  dilute  neutral  plasma  with  very  large  collision 
mean  free  time.  If  we  try  to  expand  (3.9a)  aronnd  1/r  =  0,  we  will  obtain  a  non-causal 
x{uj)  because  the  homogeneous  solution  Xh{uj)  changes  form  in  this  limiting  case.  To  see 
this,  expand  (3.8c)  aronnd  1/r  =  0  to  obtain 

(uj^Xhiuj),  (p{u))  =  (0,  (j){u)) ,  (3.15a) 
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for  all  0  G  S'.  The  form  of  the  solution  to  (3.15a)  is  given  by,  see  (A. 31), 


Xh{<^)  =  co6{uj)  +  ci6'{uj), 


(3.15b) 


which  differs  from  the  homogeneous  part  of  (3.9a)  by  a  ^^cu)  term.  Therefore,  the  full 
solution  to  =  —e^nc/m  is  given  by 

f  ^  +  co5(cu)  +  ci(5'(a;)  j  ,  (3.15c) 

where  =  ine'^ric/m  and  Up  is  called  the  angular  plasma  frequency.  To  find  the  unknown 
constants,  Cq  and  Ci,  we  use  Taylor’s  causality  requirement,  Imy  =  df  [Rey],  which  yields 
Co  =  0  and  ci  =  ivr.  Therefore,  the  causal  electric  susceptibility  is  given  by 


.X  /  1  d  ,  A 

(3.16a) 

X(t  -  t')  =  R(a,)|  =  A  “  *'>■ 

(3.16b) 

We  compare  (3.16a)  to  a  common  textbook  expression  for  the  electric  susceptibility  of  dilute 
neutral  plasma,  given  by  [16] 


oal.  1 

- -t.zrr 


(3.17) 


This  differs  from  (3.16a)  in  a  neighbourhood  of  a;  =  0.  The  difference  is  caused  by  the  fact 
that  (3.16a)  is  causal  and  valid  for  all  cc,  unlike  the  textbook  version  (3.17),  which  is  only 
valid  for  high  enough  frequencies,  namely,  cn  3>  1/r. 

The  above  example  illustrates  that  the  approximation  in  the  Fourier-domain  must  be 
done  with  care  if  we  insist  on  a  causal  electric  susceptibility.  It  is  interesting  to  note  that 
no  such  care  is  necessary  if  the  approximation  is  done  in  the  time-domain.  For  example, 
expanding  (3.9b)  around  {t  —  t')/T  =  0  immediately  yields  (3.16b),  which  is  clearly  causal. 
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3.5.1  Plasma  in  a  constant  electric  field 


Next,  let  us  consider  a  plasma  in  a  constant  electric  field  Eq  that  was  turned  on  at  t'  =  0. 
Substituting  E(t')  =  Eo0(t')  and  (3.16b)  into  (3.3a)  yields 

9 

P  Tl 

Pit)  =  V^t2Eo0(t).  (3.18) 

2m 

From  (3.16b)  and  (3.18),  we  see  that  the  electric  susceptibility  is  linear  in  time  and  that  the 
polarization  vector  is  qnadratic  in  time.  The  qnadratic  dependence  of  P(t)  on  time  signifies 
that  the  dilnte  nentral  plasma  is  accelerating  nniformly  under  the  influence  of  the  applied 
static  electric  field.  Moreover,  from  (3.3a),  we  see  that  D  grows  qnadratically  in  time  and 
that  at  t  =  0,  the  response  field  D  is  equal  to  the  applied  field  Eg. 

3.5.2  Plasma  in  a  monochromatic  electric  field 

Finally,  if  we  drive  the  plasma  with  a  monochromatic  field  of  angular  frequency  cud  that 
existed  since  t'  =  —oo,  we  will  see  another  source  of  confusion  for  stndents.  Snbstituting 
E(F)  =  Eo  cos{udt')  and  (3.16b)  into  (3.3a)  yields 

^2 

D(t)  =  eE(f)  and  e  —  1  =  Ittx  (cUd)  = - 1.  (3.19) 

Again,  if  we  replace  cud  with  u  in  (3.19)  and  ‘promote’  xi^)  to  be  a  tempered  distribution, 
then  the  ‘promoted’  xi^)  differs  from  xi^)  only  in  a  neighbonrhood  of  a;  =  0.  Also,  notice 
that  all  qnantities  in  (3.19)  are  pnrely  real.  Of  course,  D(f)  oscillates  in-phase  with  E(t) 
because  we  have  effectively  ignored  collisions  (absorption)  in  our  approximation.  If  we  let 
E  =  Egexp  (— icudt),  then  (3.19)  may  be  written  as  D(t)  =  Re(e)Re(E),  which  to  a  student 
may  confirm  an  improper  interpretation  of  (3.2). 

3.6  Damped  harmonic  oscillator 

A  damped  harmonic  oscillator  is  a  simple  model  for  the  motion  of  a  bound  electron  in  a 

dielectric.  The  eqnation  of  motion  in  the  time-domain  is  given  by 

d^r  dr  9  e  _  x 

-h  7— -h  cdgr  = - E,  (3.20a) 

dt^  dt  m 
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where  ujq  is  the  natural  angular  frequency  and  7  is  the  radiation  damping.  The  Green’s 
function  for  (3.20a)  (obtained  in  the  same  manner  and  notation  as  in  Section  3.4)  is 

(ca  —  a;+)  (ca  —  a;_)  g(a;)  =  — E(a;),  (3.20b) 

m 

where 


UJ±  =  ±s  —  1  — 
2 


and 


7 

s  =  \i 


-  .  b.,2 


(3.20c) 


If  we  assume  that  each  electron  in  the  dielectric  oscillates  with  the  same  natural  angular 
frequency,  then  using  P(ci;)  =  — enbg(ca)  yields 


{u-u+)  = 


m  ’ 


(3.21a) 


where  Uh  denotes  the  density  of  the  bounded  electrons  (assumed  to  be  constant).  From 
(3.21a),  we  see  that  x(a;)  has  two  simple  poles  in  the  complex  ca-plane.  As  we  will  soon  see, 
the  location  of  these  poles  in  the  complex  ta-plane  will  dictate  the  form  of  the  homogeneous 
solution,  Xh(i^)-  First,  consider  the  simplest  under-damped  case,  namely  when  s  >  0  and 
7  7^  0.  In  this  case,  u—u^  and  never  equal  zero  for  real  u.  Therefore,  the  homogeneous 

solution  is  simply  zero,  Xh(i^)  =  0,  and  the  full  solution  is  given  by 

1 


e^rib 


X{t  -t')  =  3^  ^  [x{uj)[  = 


m  (ca  —  a;+)  (ca  —  a;_)  ’ 
e^nb  e“2  (*“*')  sin  [s{t  —  t')] 


(3.21b) 


(3.21c) 


m  s 

Notice  that  we  didn’t  have  to  impose  Taylor’s  causality  requirement,  Imy  =  TC  [Rey],  as  it 
was  ‘automatically’  satished  by  (3.21b). 

To  get  a  better  understanding  of  (3.21c),  let’s  put  it  in  a  constant  electric  held  that  turns 
on  at  t'  =  0.  Substituting  E(F)  =  Eo0(t')  and  (3.21c)  into  (3.3a)  yields 

7  sin(st) ' 


P(i)  =  Pc 


1  —  e  2  cos{st)  + 


2s 


0(t), 


(3.22a) 


where 


e^nb 


0  — 


E, 


moo. 


0- 


(3.22b) 


0 
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From  (3.22),  we  see  that  P(t)  monotonically  increases  from  zero  at  t  =  0  to  some  maximum 
value,  and  then  oscillates  around  Pq  before  hnally  settling  at  Pq-  The  oscillations  around 
Po  remind  us  that  the  bounded  electrons  oscillate  around  the  new  equilibrium  position. 

In  the  case  of  vanishing  radiation  damping,  we  may  set  7  =  0  in  (3.21c)  to  obtain 

2 

x(t  —  t')  = - ^  sin  [a;o(t  —  t')]  Q(t  —  t'),  (3.23) 

muo 

which  is  clearly  causal.  But  if  we  set  7  =  0  in  (3.21b),  we  would  violate  causality!  To  obtain 
a  causal  7(0.;),  we  set  7  =  0  in  (3.21a)  to  obtain 


{u  -  Uo)  (uj  +  Wo)  x(ca)  = 


e^rib 

m 


(3.24a) 


From  (A. 35)  and  (A. 36),  we  see  that  the  form  of  the  full  solution  to  (3.24a)  is  given  by 

1  1 


xi(^)  = 


e^rib 

2mujn 


U  +  Uq  U  —  Uq 


+  Xhiuj) 


(3.24b) 


where 


Xh{uj)  =  boS{u  +  Wo)  +  co(5(w  -  wo).  (3.24c) 

As  in  previous  examples,  we  find  the  unknown  constants,  bo  and  Co,  via  Taylor’s  causality 
requirement,  Imy  =  TC  [Re)^,  which  yields  bo  =  —m  and  co  =  Itt.  Therefore,  the  full  casual 
solution  to  (3.24a)  is  given  by 

^  Tl\y^  1  1 

x(^)  =  7 -  — ^ - i7r5(ca  +  wq) - h  m5{u:  -  wo) 

ZtYIUJq  (jJ  UJq  CJ  —  CUq 

and  the  inverse  Fourier  transform  of  (3.24d),  of  course,  yields  (3.23).  From  the  above  exam¬ 
ple,  we  again  conclude  that  when  considering  limiting  cases  of  the  electric  susceptibility  in 
the  Fourier- domain,  we  must  be  careful  not  to  inadvertently  violate  causality.  However,  in 
the  time-domain,  we  don’t  have  to  worry  about  the  solution  not  reducing  to  a  proper  form 
in  these  limiting  cases. 


(3.24d) 
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3.7  Concluding  remarks 


In  this  paper,  we  have  taken  a  somewhat  contrarian  approach  to  the  linear  response  laws 
of  classical  electrodynamics  by  looking  at  the  response  fnnctions  in  the  time-domain.  The 
advantage  of  the  time-domain  is  that  all  qnantities  are  purely  real  and  causality  is  naturally 
expressed  in  terms  of  time.  The  disadvantage  is  that  the  response  functions  are  temporally 
non-local,  so  most  of  us  get  tired  of  writing  convolutions  on  the  blackboard  and  quickly 
slip  into  a  short-hand  mix  of  time  and  frequency/ Fourier  domain  notations  that  can  confuse 
students  profoundly. 

While  it  is  perfectly  reasonable  to  avoid  complications,  such  as  dispersion,  in  introductory 
physics  courses,  by  the  time  students  are  in  their  third  or  fourth  year  of  physics  study,  it  is 
important  to  expose  them  to  the  fundamental  principles  associated  with  a  classical,  macro¬ 
scopic  picture  of  matter.  In  particular,  we  believe  that  the  following  should  be  emphasized: 

•  Causality  is  easy  to  enforce  in  the  time-domain,  but  the  constitutive  relations  are 
non-local  in  time  and  involve  convolution  integrals. 

•  The  constitutive  relations  are  mathematically  simple  in  the  Fourier-domain,  but  causal¬ 
ity  is  given  by  the  Kramers-Kronig  relations  (the  Hilbert  transform  pair). 

•  We  should  make  it  clear  whether  we  are  really  transforming  into  the  Fourier-domain, 
or  whether  we  are  assuming  a  monochromatic  source  in  a  time-domain  experiment. 
Very  often  the  results  look  the  same,  but  as  we  have  shown,  confusing  the  two  can  lead 
to  serious  misunderstandings. 

•  Finally,  it  should  be  emphasized  that  Maxwell’s  equations  (in  the  time-domain)  are 
purely  real  and  involve  only  purely  real  quantities.  The  Fourier  transformation  pro¬ 
motes  variables  to  the  complex  plane.  As  teachers,  we  should  be  careful  when  speaking 
of  the  real  and  imaginary  parts  of  the  response  function,  and  preface  our  remarks  with 
a  note  that  we  are  working  in  the  non-physical,  but  highly  useful,  Fourier-domain. 
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A.l  Appendix  A  -  Distribution  theory 

The  ordinary  (classical)  functions  may  be  thought  of  as  a  mapping  between  two  sets 
of  numbers.  We  can  extend  the  idea  of  an  ordinary  function  by  considering  a  mapping 
(functional,  if  you  will)  between  a  set  of  functions  and  a  set  of  numbers.  The  physical 
reason  for  extending  the  idea  of  an  ordinary  function  lies  in  our  inability  to  experimentally 
measure  a  function  at  a  point,  e.g.,  see  [21,  pp.  1-2],  [22].  Let  f{x)  represent  temperature 
at  some  point  x.  To  measure  the  temperature  at  that  point,  we  place  a  thermometer  there 
to  obtain  a  value  for  f{x),  but  do  we  actually  obtain  a  value  of  the  temperature  at  point 
x7  The  bulb  of  the  thermometer  has  hnite  size;  thus,  what  we  measure  is  an  average 
temperature  around  the  point  x.  Mathematically  an  average  is  a  weighted  sum,  and  our 
measurement  only  reveals  the  value  Ti  =  ff(x)</)i(x)  dx,  where  <fi(x)  is  essentially  zero  away 
from  the  bulb  of  the  thermometer.  If  we  make  another  measurement  of  the  temperature 
at  point  X  using  a  different  thermometer,  then  we  would  measure  T2  =  ff(x)<f2(x)dx, 
and  hopefully,  the  ‘true’  temperature  T  ^  (Ti  +  T2)  /2.  The  above  discussion  is  meant  to 
serve  as  a  physical  motivation  for  defining  what  we  will  call  generalized  functions  as  certain 
linear  functionals.  Generalized  functions  are  also  called  distributions,  and  we  will  use  both 
terms  interchangeably.  Our  presentation  of  the  generalized  functions  closely  follows  that 
of  Strichartz  [21].  A  reader  interested  in  a  more  detailed  study  may  also  hnd  [23-26]  helpful. 

For  our  purposes,  it  will  be  sufficient  to  consider  only  a  special  class  of  distributions, 
namely,  the  tempered  distributions.  Before  we  can  formally  dehne  tempered  distributions, 
we  must  hrst  dehne  a  set  of  ‘good’  functions.  This  set  of  ‘good’  functions  is  also  called 
the  space  of  test  functions.  For  the  tempered  distributions,  the  space  of  test  functions, 
denoted  by  S,  contains  all  real  or  complex-valued  functions  0(t)  that  are  classically  inhnitely 
differentiable  and,  along  with  all  its  derivatives,  vanish  at  inhnity  faster  than  the  reciprocal 
of  any  polynomial.  For  example,  any  function  of  the  form  '^n=o  ^^t"' exp  {—t"^)  belongs  to 
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S.  We  are  now  ready  to  define  the  class  of  tempered  distributions,  denoted  by  S',  as  all 
continuous^  linear  functionals  on  S.  A  linear  functional  /  on  S  is  a  rule  by  which  we  assign 
to  every  test  function  (fit)  a  (real  or  complex)  number  denoted  by  such  that  the 

identity  (/,  cicfi  +  02^2)  =  Ci  (/,  ^i)  +  C2  (/,  ^2)  is  satished  for  arbitrary  test  functions  4>i 
and  02  and  (real  or  complex)  numbers  ci  and  C2.  The  terminology  and  notation  used  for 
distributions  can  be  confusing  at  times  because  the  phrase  ‘function  /  or  even  generalized 
function  (distribution)  /’  may  refer  to  /  itself  or  to  the  value  of  (/,  0).  In  other  words,  no 
distinction  is  made  between  a  distribution  and  a  ‘function’  from  which  the  distribution  was 
obtained.®  To  make  the  notion  of  tempered  distributions  more  concrete  let’s  consider  a  few 
simple  examples.  Let’s  hnd  a  distribution  dehned  by  /  =  0'(t),  i.e. 


7f  = 


d0(t) 

dt 


(j){t)  dt  = 


^  ^  dt 


d0(t) 

dt 


dt  =  -  (0(oo)  -  0(0))  =  0(0), 


(A.l) 

(A.2) 


where  we  integrated  by  parts  in  (A.l)  and  the  integrated  terms  vanished  because  0  is  a 
‘good’  function,  i.e.  0  G  S'.  By  comparing  (A.2)  to  the  sifting  property  of  the  Dirac  delta 
function 


/OO 

(5(t)0(t)dt  =  0(0),  (A.3) 

•OO 

we  conclude  that  the  (generalized)  derivative  of  the  Heaviside  step  function  equals  the  Dirac 
delta  function.  We  can  even  differentiate  (in  a  distributional  sense,  of  course)  more  compli¬ 
cated  functions.  Let  /  =  g{t)5'{t),  where  g{t)  is  a  polynomial  of  any  degree;  then 


^The  definition  of  continuity  of  linear  functionals  is  rather  technical  and  not  necessary  for  our  purposes. 
^Strictly  speaking,  this  is  an  abuse  of  terminology  but  it  is  so  common  that  one  must  be  aware  of  it.  Moreover, 
one  often  speaks  of  generalized  functions  (distributions)  as  if  they  were  proper  functions,  e.g.,  the  Dirac 
delta  function  5{t). 
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Integrating  by  parts  and  noticing  that  the  prodnct  of  g(t)  and  0(t)  is  still  in  S  (so  that  the 
integrated  terms  vanish),  yields 

+ ^^<‘0  * 

=  _^(O)0'(O)- ^'(0)0(0).  (A.5) 

We  can  write  (A.5)  in  a  more  standard  form  that  doesn’t  involve  the  derivatives  of  (p{t). 
Using  (A. 3)  and  noting  that 

we  obtain 

-  {g'{0)6{t),(j){t)) .  (A.6) 

It  is  a  very  common  abuse  of  notation  to  ‘drop’  the  ( ,  )  brackets,  along  with  4>{t),  and  write 

(A.6)  simply  as 

g{t)5'{t)=g{0)5'{t)-g'{0)5{t).  (A.7) 

In  particular,  if  we  let  g(t)  =  t  in  (A. 4),  then  (A.7)  yields  tS'(t)  =  — 5(t);  not  just  zero  as 
one  might  have  naively  expected.  An  alert  reader  may  have  noticed  that  in  the  derivation 
of  (A.7),  we  never  used  the  assumption  that  g{t)  is  a  polynomial;  all  that  the  derivation 
required  was  g(t)(l)(t)  G  S.  While  it  is  dehnitely  true  that  g{t)(l){t)  G  S  when  g{t)  is  a 
polynomial,  requiring  g{t)  to  be  a  polynomial  is  an  unnecessary  restriction.  In  other  words, 
(A.7)  holds  for  any  function  g{t)  as  long  as  g(t)(l)(t)  G  S.  For  example,  g(t)  could  be  sin(t), 
but  it  cannot  be  exp  (t^)  because  then  g{t)(l){t)  ^  S  and  the  integrated  terms  will  not  vanish. 
We  considered  this  example  in  such  detail  because  we  will  have  numerous  opportunities  to 
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use  (A. 7)  in  appendix  A. 1.2. 

A.  1.1  The  Fourier  transform  of  tempered  distributions 

The  key  idea  in  generalizing  the  notation  of  a  derivative  is  to  move  the  derivative  from  a 
function  (generalized  function)  onto  a  set  of  ‘good’  functions,  namely,  0(t)  G  S.  Moreover,  we 
saw  that  by  integrating  by  parts  enough  times,  every  (generalized)  function  had  a  derivative, 
because  the  space  S  is  composed  of  classically  inhnitely  differentiable  functions.  We  will 
use  this  ‘moving  idea’  (adjoint  operator)  to  define  the  Fourier  transform  of  a  tempered 
distribution  f(t).  By  the  Fourier  transform  of  f(t)  G  S',  denoted  by  f{u)  or  by  T[/(t)],  we 
mean 

(3"  [fit)] ,  =  {fit),  S'  [0(w)]) ,  (A.8) 

where 

/oo  /  roo  \ 

fit)  I  0(a;)e+'‘^*  dw  J  dt.  (A.9) 

OO  \j  —OO  / 

By  the  inverse  Fourier  transform  of  /(ca),  denoted  by  fi^)  ,  we  mean 

fiuj)  ,0(t)^  =  ^/(a;),T“^  fiit)  (A.IO) 

where 

(7^,3""^  fit))  =  j  7^1^^^  fiit)e-'‘^^d?j  du.  (A.ll) 

By  changing  the  order  of  integration  in  (A.9)  and  (A.ll),  we  see  that  (A.9)  and  (A.ll) 
are  indeed  consistent  with  the  classical  Fourier  transform  pair.  The  underlining  reason  for 
dehning  the  Fourier  transform  pair  by  (A.8)  and  (A.IO)  lies  in  the  fact  that  if  0  G  A  then 
0  is  also  in  A;  the  converse  is  also  true.  Moreover,  it  can  be  shown  that  /  G  S"  if  and  only 
if  /  G  S'.  For  a  proof  of  these  and  related  matters,  see  [21,  chapter  3],  [24,  chapter  6],  [25, 
chapter  7]  and  [26,  chapter  2].  Now,  we  will  make  (A.8)-(A.ll)  more  concrete  by  considering 
a  few  simple  examples. 
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For  our  first  example,  we  will  compute  the  generalized  Fourier  transform  of  S{t),  i.e. 
(3'[(5(t)]  ,0(a;)).  Substituting  S(t)  into  (A. 9)  and  changing  the  order  of  integration  yields 


0(cn) 


dt  den  =  /  0(a;)  den  =  (1,  0(en)) . 


(A.12) 


Thus,  we  see  that  T[(5(f)]  =  1.  Moreover,  by  taking  the  inverse  Fourier  transform,  we  obtain 
the  famous  integral  representation  of  the  Dirac  delta  function,  namely. 


den. 


(A.13) 

It’s  worth  stressing  that  (A.13)  should  be  interpreted  in  a  distribution  sense  and,  strictly 
speaking,  writing  (A.13)  as  we  did  is  an  abuse  of  notation.  However,  such  abuses  of  notation 
are  very  common  in  physics;  e.g.,  see  the  famous  graduate  electrodynamics  textbook  [27]. 

As  another  simple  example,  consider  the  generalized  inverse  Fourier  transform  of 
2n6{u  —  Uq),  where  cno  is  a  real  number;  i.e.  27i6{u  —  Uq)  when  cno  G  M. 

Substituting  27i6{u  —  cno)  into  (A. 11),  then  changing  the  order  of  integration  and  integrating 
over  cn  yields 


0(t)e 


dt  =  ( 


(A.  14) 


Thus,  we  see  that  (5(cn  —  cno)  =  exp(— icnot). 

For  our  last  example,  let  us  compute  the  Fourier  transform  of  the  nth  generalized  deriva¬ 
tive  of  f{t)  G  S' .  From  (A. 8),  we  have 


T 


,0(cn)  )  = 


r-oo  /  jn  \  ^ 


(A.15) 


Performing  integration  by  parts  n-times  on  the  right-hand  side  (where  the  integrated  terms 
vanished  because  0  G  S'),  yields 


(-1)”  /  f{t)  (  J  dt  =  (-1)”  J  /(t)T  [(Ticn)"^  0(cn)]  dt 

poo 

=  /  3'[/(^)]  (-it^)''0(^^)dcn. 


(A.16) 
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where  we  obtained  (A. 16)  by  using  the  dehnition  (A. 8).  Finally,  comparing  (A. 16)  with  the 
left-hand  side  of  (A.  15)  yields 


5 


(-kj)”  ?!/(«)] . 


(A.17) 


Notice  that  in  this  example,  we  have  used  the  dehnition  (A. 8)  twice]  Strichartz  [21,  pp.  49-50] 
appropriately  refers  to  this  as  ‘dehnition  chasing’. 

No  discussion  of  the  Fourier  transform  of  tempered  distributions  would  be  complete 
without  the  Fourier  transform  of  a  convolution  integral.  For  our  purposes,  it  will  be  sufficient 
to  only  consider  convolution  of  the  tempered  distribution  /  with  a  hxed  element  ip  from  the 
space  of  ‘good’  functions.  Let  f  E  S'  and  ip  E  S]  then,  by  the  convolution  of  /  with  ip, 
denoted  by  (/  *  ip)  (t),  we  mean 


(f*^P)(t)=  f{t-t')ip{t')dt'.  (A.18) 

J  —  OO 

By  a  simple  change  of  variables,  we  see  that  convolution  is  commutative,  i.e.  (/  *  ip)  {t)  = 
{ip  *  /)  (t).  Convolution  dehnes  an  inhnitely  diherentiable  function  and  thus  can  be  viewed 
as  a  ‘smoothing’  process.  To  see  this,  let  h{t)  =  {f  *ip)  {t)  then 


=  ^  [(^  *  /)  W]  =  j  (^-20) 

From  (A. 20),  we  conclude  that  all  derivatives  of  h{t)  exist  in  the  classical  sense  because 
Ip  E  S,  and  from  (A.  19)  we  see  that  it  doesn’t  matter  how  ‘rough’  (e.g.,  6'{t)  is  ‘rougher’ 
than  S{t))  the  distribution  is.  In  passing,  we  note  another  useful  property  of  (A.18),  namely 
that  its  Fourier  transform  corresponds  to  multiplication  in  the  Fourier-domain,  i.e. 


(A.21) 


A.  1.2  Support  and  structure  of  tempered  distributions 

We  often  want  to  speak  about  the  local  properties  of  distributions  as  if  the  distributions 
were  ordinary  (classical)  functions.  When  speaking  about  an  ordinary  function  f{t),  the 
statement  ‘/(^)  has  a  value  of  f{ti)  when  t  =  ti  has  meaning,  but  the  statement  is  nonsense 
if  f{t)  is  a  distribution.®  However,  we  can  identify  a  set  of  points  where  distribution  /  is  non¬ 
zero.  Loosely  speaking,  this  set  of  points  is  known  as  the  support  of  /.  The  formal  dehnition 
of  support  requires  us  to  dehne  where  a  distribution  /  is  zero.  We  say  a  distribution  f(t)  is 
zero,  fit)  =  0,  on  an  open  interval  (a,  6)  if  {fit), 'fit))  =  0  for  every  inhnitely  differentiable 
test  function  if  it)  that  vanishes  in  a  neighbourhood  of  every  point  not  in  (a,  6)  interval.  For 
example,  if  fit)  =  0  on  (—1, 1),  then  iff)  vanishes  outside  the  (— l-|-e,  1  — e)  interval  for  some 
e  >  0.  We  now  formally  dehne  the  support  of  a  distribution  ff),  denoted  by  supp  [/(t)], 
as  the  complement  of  the  set  of  points  t  such  that  ff)  =  0  in  a  neighbourhood  of  t.  For 
example,  if  ff)  =  0  on  (— cxd,cx)),  then  supp  [/(t)]  is  the  empty  set  and,  as  a  less  trivial 
example,  supp  [5(t)]  =  {0}.  Moreover,  we  will  show  that 

supp  ={0}-  (A.22) 

Consider  any  open  interval  Ji  that  does  not  contain  the  point  t  =  0;  then, 

because  iff)  vanishes  in  a  neighbourhood  of  t  =  0.  Of  course,  for  any  open  interval  I2  that 
does  contain  the  point  t  =  0,  the  derivatives  of  if  don’t  vanish  at  the  point  t  =  0  for  every 
test  function  if.  Thus,  we  see  that  all  derivatives  of  Sf)  have  the  same  point-support. 

The  above  discussion  was  necessary  to  understand  the  following  ‘structure’  theorem.  A 
tempered  distribution  ff)  with  supp  [ff)]  =  {to}  must  be  of  the  form  [21,  pp.  82-88] 

/w  =  (5(t-to),  (A.23) 

n=0 

^Recall  that  we  have  defined  distributions  only  by  their  action  on  the  space  of  test  functions. 
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where  the  coefficients  Cn=o,...,N  are  complex  numbers,  i.e.  G  C  for  n  =  0, . . . ,  iV.  In 
other  words,  any  tempered  distribution  with  a  point-support  may  be  expressed  as  a  finite 
linear  combination  of  the  Dirac  function  and  its  derivatives;  this  is  a  powerful  statement! 
In  the  next  paragraph,  we  will  show  how  we  can  use  (A. 23)  to  solve  ‘algebraic’  equations  in 
a  distributional  sense.  These  types  of  equations  are  frequently  encountered  when  we  solve 
differential  equations  by  the  Fourier  transform  technique. 

As  our  first  example,  consider  the  following  equation, 

((f  -  fo)  f{t),  =  (1,  ,  (A.24) 

for  all  (j){t)  G  S.  Before  we  hud  the  unknown  tempered  distribution  f{t),  we  note  that  it  is 
customary  to  abuse  the  notation  and  write  (A.24)  simply  as 

{t-to)f{t)  =  l.  (A.25) 

Naively,  we  might  expect  that 

/p(«)  =  (A.26) 

would  solve  (A.25),  but  this  is  only  the  particular  part  of  the  solution.  We  could  have  a 
tempered  distribution  /h(t)  with  supp  [fh(t)]  =  {to}  such  that  (t— to)/h(t)  =  0,  and  therefore, 
fit)  =  /p(t)  -|-/h(t)  would  also  satisfy  (A.25).  We  refer  to  /h(t)  as  the  homogeneous  solution 
and,  in  light  of  the  structure  theorem  in  the  previous  paragraph,  we  know  it  must  be  of  the 
form 

hit)  =  ^^€^^6  it -to).  (A.27) 

ri=0 

Substituting  (A.27)  into  ((t  —  to)/h(t), 0(t))  =  (O,0(t))  and  integrating  by  parts  until  the 
derivatives  only  appear  on  yields 

-Clfiito)  +  2c20'(to)  -  3C30"(to)  H - h  i-1)^ NcN(j)^^~^\to)  =  0. 
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The  above  equation  must  hold  for  all  (j)  E  S.  Thus,  the  coefficients  must  vanish  independently, 
i.e.  Cn  =  0  for  n  =  1,2, . . . ,  N .  Therefore,  the  homogeneous  solution  is  given  by 


fh{t)  =  CoS{t -to),  (A.28) 

and  the  full  solution  to  (A. 25)  (or  more  formally,  to  (A. 24))  is  given  by 

f(t)  =  - — —  +  Co6{t  —  to).  (A. 29) 

t  —  to 

Loosely  speaking,  from  (A. 29)  we  see  that  we  can  divide  by  zero,  provided  we  add  an  ap¬ 
propriate  tempered  distribution  with  a  point-support.  In  general,  using  the  same  procedure 
as  above,  we  can  show  that  the  distributional  solution  to  {t  —  to)^  f{t)  =  1  is  given  by 
fit)  =  fpit)  +  fhit),  where 

/p(t)  =  (^-30) 

fh{t)  =  CoS{t  —  to)  +  CiS\t  —  to)  +  ■  ■  ■  +  ^\t  —  to).  (A. 31) 


For  our  second  and  last  example,  consider  (in  a  distributional  sense,  of  course)  the  fol¬ 
lowing  equation, 

{t-h){t-t2)f{t)  =  l,  (A.32) 

where  ti  7^  ^2-  From  our  previous  example,  we  expect  the  solution  to  be  of  the  form 

N  M 

^ <*-*■>  +  E ■  (A'®) 

'  n=0  m=0 

Substituting  (A. 33)  into  (A.32)  (of  course,  we  actually  mean  {{t  —  ti)  (t  —  t2)  f  (t) ,  (pit))  = 

(1,  0(t)),  for  all  4)  E  S)  and  integrating  by  parts  until  the  derivatives  only  appear  on  0  yields 
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From  (A. 34),  we  see  that  =  0  for  n  =  1, . . . ,  A^  and  =  0  for  m  =  1, . . . ,  M.  Therefore, 
the  full  solution  to  (A. 32)  is  given  by  f(t)  =  fp(t)  +  fh(t),  where 


'^‘’***  ((-(,)((  -(2)' 

fh{t)  =  bo6  {t  -  ti)  +  cqS  {t  -  t2) . 


(A.35) 

(A.36) 


Notice  that  if  ti  =  t2  then  (A.36)  does  not  yield  the  correct  solution,  which  is  given  by  (A. 31) 
(with  n  =  2  and  ti  =  t2  — t  to).  The  reason  for  this  ‘discrepancy’  is  because  our  conclusion 
from  (A. 34),  namely,  that  bn=i,...,N  =  0  and  Cm=i,...,M  =  0,  is  not  valid  if  ti  =  t2.  In  other 
words,  we  must  be  very  careful  when  dealing  with  distributional  solutions  in  limiting  cases 
such  as  t2  — )■  ti-  In  the  body  of  the  paper,  these  limiting  cases  arise  when  we  consider 
vanishing  absorption. 
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CHAPTER  4 


MEASURING  THE  VOID:  THEORETICAL  STUDY  OF  SCATTERING  BY  A 

CYLINDRICAL  ANNULUS 

A  paper  published  in  the  Journal  of  Quantitative  Spectroscopy  &  Radiative  Transfer. 

Alex  J.  Yuffa*,  John  A.  Scales 

Department  of  Physics,  Colorado  School  of  Mines,  Golden,  CO  8O4OI,  USA 
*Primary  researcher  and  author.  E-mail:  ayuffa@gmail.com 

4.1  Abstract 

In  this  paper,  we  analyze  a  monochromatic  plane  wave  scattering  from  an  inhnite  ho¬ 
mogeneous  cylindrical  annulus.  In  particular,  we  study  the  effect  that  the  inner  part  of  the 
cylindrical  annulus  (cylindrical  void,  if  you  will)  has  on  the  scattered  held.  This  is  done  by 
isolating  the  cylindrical  void’s  contribution  to  the  scattered  held.  We  show  that  if  the  cylin¬ 
drical  void  is  small,  then  its  contribution  to  the  scattered  held  may  be  approximated  by  the 
“screened  cylindrical  void”  (SCV)  approximation.  We  hrst  develop  the  SCV  approximation 
in  a  physically  intuitive  manner,  and  then  show  that  it  could  also  be  obtained  in  a  more 
mathematically  rigorous  manner.  Numerical  results  comparing  the  SCV  approximation  to 
the  exact  solution  are  also  presented. 

4.2  Introduction 

Consider  a  monochromatic  plane  wave  scattering  from  an  inhnitely  long  homogeneous 
and  isotropic  cylindrical  annulus  with  outer  radius  ri  and  inner  radius  r2,  see  Figure  4.1(a). 
Let  Cl  denote  the  permittivity  of  the  space  surrounding  the  cylindrical  annulus  and  let  62 
denote  the  permittivity  of  the  cylindrical  annulus  itself,  r2  <  r  <  ri.  Let  us  refer  to  the 
region  of  space  inside  the  cylindrical  annulus  as  the  “cylindrical  void”  and  ask  what  effect 
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the  cylindrical  void  has  on  the  scattered  £eld(s)  outside  the  cylindrical  annulus.  If  one  were 
to  experimentally  investigate  this,  one  would  do  the  following: 


(a)  measure  the  total  field  V^^\r^9)  outside  the  cylindrical  annulus  (r  >  ri); 

(b)  measure  the  total  field  6)  outside  an  identical  “host  cylinder;”  i.e.,  a  cylinder  of 

radius  ri  and  permittivity  62,  as  illustrated  in  Figure  4.1(b); 

(c)  compute  the  difference  between  the  two  fields  in  (a)  and  (b): 

l^(sca)(^^^)  _  v^^\r,e)-U^^\r,e).  (4.1) 


Following  the  above  procedure,  contains  the  effect  that  the  cylindrical  void 

had  on  the  scattered  held.  In  this  paper,  we  show  that  can  be  approximated 

by  the  scattered  held  produced  by  the  cylindrical  void  when  a  plane  wave  from  a  region  of 
space  with  a  permittivity  of  62  is  incident  on  it.  This  approximation  holds  if  the  “screening 
effect”  (discussed  in  Section  4.3)  of  the  cylindrical  annulus  is  properly  accounted  for,  and 
if  the  cylindrical  void  is  sufficiently  small.  We  refer  to  this  approximation  as  the  sereened 
cylindrical  void  (SCV)  approximation.  Furthermore,  we  investigate  the  rate,  denoted  by 
at  which  the  energy  is  extinguished  (depleted)  by  the  cylindrical  void  from  the  total 
held  outside,  U^^\r,9),  the  host  cylinder. 

To  the  best  of  our  knowledge,  the  SCV  approximation  and  its  physical  interpretation 
(see  Section  4.3)  has  not  been  previously  considered  in  the  literature.  In  order  to  make 
the  paper  accessible  to  the  widest  possible  scientific  community,  we  use  the  well-known 
Lorenz-Mie  theory  [1-4]  to  derive  the  SCV  approximation.  However,  we  do  note  that  our 
intuitive  derivation  of  the  SCV  approximation,  which  is  presented  in  Section  4.3,  is  physically 
guided  by  the  Debye  series  expansion  [5].  In  short,  the  Debye  series  expansion  consists  of 
re-expressing  each  Mie  scattering  coefficient  in  terms  of  an  infinite  series  called  the  Debye 
series.  Each  term  in  the  Debye  series  may  be  physically  interpreted  in  terms  of  the  number 
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Figure  4.1:  The  cross-sectional  view  of  the  cylindrical  scattering  objects  is  shown.  The  origin 
of  the  coordinate  system  (r,  6),  where  — tt  <  6*  <  tt,  is  concentric  with  the  cylindrical  objects. 
In  each  panel,  the  region  is  denoted  by  a  boxed  number  and  the  permittivity  of  each  region 
is  also  indicated.  For  example,  region  three,  r  <  r2,  in  panel  (a)  has  a  permittivity  of  ei  and 
region  one,  r  >  r2,  in  panel  (c)  has  a  permittivity  of  62. 


of  reverberations  the  wave  has  experienced.  A  reader  interested  in  the  use  of  the  Debye 
series  expansion  in  the  context  related  to  this  paper,  namely,  plane  wave  scattering  by  a 
multilayered  cylinder,  may  consult  [6,  7]  and  references  therein. 

Although  we  do  not  explicitly  consider  many  diverse  areas  of  science  where  the  scattering 
by  a  cylindrical  void  is  important  (e.g.,  see  [3,  4]),  we  would  like  to  mention  one,  namely, 
localization.  Fifty  years  after  the  pnblication  of  Anderson’s  seminal  work  [8],  localization 
continues  to  be  a  thriving  area  of  research  [9]  in  theoretical  and  experimental  physics.  Local¬ 
ization  of  millimeter /snbmillimeter  electromagnetic  waves  is  particularly  interesting  because 
both  the  amplitude  and  the  phase  of  the  electromagnetic  held  can  be  easily  measured  with  a 
vector  network  analyzer  [10].  At  these  wavelengths,  the  preparation  of  disordered  samples  is 
also  inexpensive  and  straightforward  with  standard  compnter-nnmerically-controlled  (CNC) 
milling  techniqnes.  A  sample  may  be  prepared  by  drilling  small  holes  in  a  large  Tehon  (nltra 
low-loss  material)  cylinder.  Further,  by  illuminating  the  sample  from  the  side  and  pntting  it 
on  a  rotational  stage,  we  can  generate  essentially  arbitrary  realizations  of  the  same  random 
disorder.  When  the  number  of  small  scatterers  is  large,  say,  over  1000,  then  what  is  impor- 
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tant  is  the  rate  at  which  the  scatterer  extinguishes  the  energy  from  the  incident  field,  rather 
than  the  geometrical  shape/size  of  each  individual  scatterer  [11,  12],  Thus,  the  physical 
insight  into  scattering  by  a  single  small  cylindrical  void  discussed  in  this  paper  may  be  of 
beneht  in  understanding  the  experimental  model  described  above. 

Throughout  this  paper,  we  will  use  the  Gaussian  unit  system,  and  we  will  assume  that  all 
helds  are  harmonic  in  time  with  a  exp(— icat)  time  factor,  where  u  is  the  angular  frequency. 
Furthermore,  we  will  assume  that  all  helds  are  polarized  in  the  positive  z-direction.  The 
positive  z-direction  is  out  of  the  page  in  Figure  4.1.  All  media  considered  in  this  paper  are 
assumed  to  be  non-magnetic,  and  ei  is  assumed  to  be  purely  real. 

4.3  Intuitive  derivation  of  the  SCV  approximation 

In  this  section,  a  physically  intuitive  derivation  of  the  SCV  approximation  is  presented. 
The  derivation  is  organized  as  follows.  First,  we  imagine  a  unit  plane  wave  inci¬ 

dent  from  region  one  onto  the  cylindrical  void  shown  in  Figure  4.1(c).  Then,  we  compute 
the  scattered  held  in  region  one  generated  by  the  scattering  of  u^^^^\r,9)  from 

the  cylindrical  void.  Second,  to  account  for  the  screening  ehect  of  the  cylindrical  annu¬ 
lus,  we  use  the  previously  found  scattered  held  as  the  incident  (primary)  held, 

i.e.,  w^^^^\r,9)  =  u^^^^\r,9),  originating  from  the  center  of  the  host  cylinder  shown  in  Fig¬ 
ure  4.1(b).  Finally,  we  compute  the  total  held  9)  in  region  one  shown  in  Figure  4.1(b) 

and  physically  interpret  the  terms  contained  in  it  to  obtain  an  approximation  to  9), 

see  (4.1). 

Let  us  note  that  all  helds  in  this  paper  satisfy  the  two-dimensional  (2D)  Helmholtz  equa¬ 
tion.  The  radial  solution  of  the  2D  Helmholtz  equation  is  composed  of  a  linear  combination 
of  integer  order  Bessel  functions  of  the  hrst  and  second  kind,  which  we  denote  by  Jn{0 
and  Yn{^),  respectively.  The  Bessel  functions  Jn{0  Yn{^)  also  satisfy  the  Wronskian 
relationship  [13],  namely, 

J„+l(0V^(0  -  ^n(0>^n+l(0  =  4-  (4-2a) 
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Also,  JniO^  ^n(0  and  the  Hankel  function  of  the  first  kind,  which  we  denote  by  Hn{^)  = 
JniO  +  satisfy  the  recurrence  relation  [13], 

H  71 

=  (4.2b) 

where  'h  denotes  J,  Y  or  H .  Lastly,  we  note  the  Jacobi-Anger  expansion  of  a  plane  wave 
[13],  namely, 

OO 

ggcose  ^  '^gnYJn{C)  cos  ine) ,  (4.2c) 

where  gn  denotes  the  Neumann  factor:  go  =  1  and  gn  =  2  for  n  >  1. 

Returning  to  the  scattering  of  the  unit  plane  wave  from  the  cylindrical  void  shown  in 
Figure  4.1(c),  let  the  incident  wave  be  6)  =  exp(iA;2^"  cos 6*),  where  /c2  =  is  the 

wavenumber  and  c  is  the  speed  of  light  in  a  vacuum.  Then,  the  field  in  region  two 
and  the  total  held  in  region  one  decomposed  as  u^^\r,9)  =  +  u^^^^\r,9),  may  be 

written  as 


-^(inc)(^,^)- 

OO 

Jn{k2r) 

^(sca)j^,  Q-j 

SnHn{k2r) 

1 

n=0 

Jn-Jnikir)  _ 

cos(n6'). 


(4.3) 


where  ki  =  ^/eiu/c.  In  writing  (4.3),  we  used  the  Jacobi-Anger  expansion  (4.2c)  to  rewrite 
exp(iA;2’"  cos  0)  as  an  inhnite  sum,  imposed  the  Sommerfeld  radiation  {outgoing  cylindrical 
wave)  condition  on  9),  and  required  9)  to  be  regular  (hnite)  at  r  =  0.  To  hnd 

the  unknown  coefficients  in  (4.3),  we  require  that  the  electric  held  and  its  normal  derivative 
be  continuous  across  the  r  =  r2  interface,  i.e.. 


.y(l)  ^  .^(2) 


and  on  r  =  r2. 


(4.4) 


dr  dr 

to  obtain  a  system  of  linear  equations.  Solving  this  system  of  linear  equations  for  Sn  and 
using  (4.2b)  to  simplify  the  result,  yields 

Jn+i{kir2)Jn{k2r2)  -  KJn{kir2)Jn+i{k2r2) 


Sn  = 


Jn+l{kir2)Hn{k2r2)  -  KJn{kir2)Hn+l{k2r2)' 


(4.5a) 
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where  k  =  k2/ki,  n  G  Z"*",  and  Z"*"  denotes  the  set  of  all  nonnegative  integers.  It  is  convenient 
to  introdnce  cnrly  bracket  notation,  $(^7)},  by  which  we  mean 


{vl/,+l(0;  ®n(h)}  =  (d^n+l(0<^>n(h)  “  (O^n+1  (h) )  • 


For  example,  (4.5a)  in  the  cnrly  bracket  notation  reads  as 

.  _  {Jn+l(fcir2);Jn(fc2r2)}  +  .  . 

{Jn^.{k,r2YH^{k.r2)y 

Having  found  the  expansion  coefficients  of  the  scattered  wave  6^),  we  are  now  ready  to 


(4.5b) 


see  how  they  should  be  modihed  in  order  to  account  for  the  screening  effect  of  the  cylindrical 
annulus. 

Imagine  a  “line-source”  embedded  in  the  center  of  the  host  cylinder  shown  in  Fig¬ 
ure  4.1(b).  We  take  the  held  produced  by  the  line-source  to  be  equal  to  in 

(4.3).  If  we  use  this  held  as  the  incident  held,  i.e.,  0),  then  the  to¬ 

tal  held  w^‘^\r,9)  inside  the  host  cylinder  (region  two  in  Figure  4.1(b))  may  be  written  as 
9)  =  9)  +  9),  where 


y^(sca)  9)  =  gj^  ldnJn{.k2r)  COS  {n9) . 


(4.6a) 


Notice  that  in  (4.6a),  we  required  w^^^^\r,9)  to  be  regular  at  r  =  0.  This  requirement  is 
necessary  because  we  are  essentially  treating  the  cylindrical  void  as  a  line-source  in  this 
paragraph.  The  held  outside  the  host  cylinder  (region  one  in  Figure  4.1(b)),  6*),  must 


satisfy  the  Sommerfeld  radiation  condition  and  thus,  it  is  given  by 


w^^\r,9)  =  y^  gni^anHnikir)  cos  (7^6^) . 


(4.6b) 


Imposing  the  boundary  conditions  and  {d/dr)  =  {d/dr)  on  r  =  ri,  then 

solving  the  resultant  linear  system  for  and  using  (4.2a)  with  (4.2b)  to  simplify  the  result 


yields 


Tf/ciri  {Hn+i{kiri)]  Jn{k2ri)} 


Sn,  n  G  Z"*". 


We  physically  interpret  the  term  in  parentheses  in  (4.7)  as  the  screening  effect  of  the  cylin¬ 
drical  annulus  on  the  scattered  wave  generated  by  the  cylindrical  void.  The  coefficients 
are  not  quite  the  correct  ones  to  use  in  (r,  9)  because  they  do  not  contain  the  screening 

effect  that  the  cylindrical  annulus  had  on  the  incident  wave.  A  moment’s  thought  reveals 
that  this  screening  effect  had  to  be  the  same  as  the  screening  effect  on  the  scattered  wave. 
Thus,  the  expansion  coefficients  should  be  given  by  (4.7)  with  the  parenthesis 

term  squared.  Therefore,  is  approximately  given  by 


ri=0 


-2i 


nkiTi  {Hn+i{kiri)]  Jn{k2ri)] 


2 

5nHn{kir)  cos  {n6) , 


where  the  5n  coefficients  are  given  by  (4.5). 


(4.8) 


4.4  Rigorous  derivation  of  the  SCV  approximation 


In  this  section,  we  present  a  rigorous  derivation  of  by  directly  computing 

U^^\r,9)  and  V^^\r,9)  (recall  the  bullet  list  of  Section  4.2).  Once  the  exact 
is  found,  we  show  that  it  is  approximately  equal  to  (4.8)  if  kir2  1  and  \k2\r2  1. 

Furthermore,  a  numerical  illustration  of  the  SCV  approximation  is  also  presented. 

If  a  plane  wave,  =  exp(iA;ir  cos  0),  is  incident  on  the  host  cylinder  shown  in 

Figure  4.1(b),  then  by  proceeding  as  in  paragraph  three  of  Section  4.3,  the  total  held  in 
region  one  is  9)  =  7/*'“'’’^(r,  9)  +  9),  where  the  scattered  held  is 

00 

U^^^^\r,9)  =  gni"' Hnjkir)  cos{n9)  (4.9a) 

n=0 

with 


.(he)  _  {Jn+i{kiri);  Jn{k2ri)}  ^  ^ 

r  TT  /  7.  „  \ .  T  /  7.  „  \  I  ’  Tl  ^  /L  , 


{Hn+i{kiri)-,  Jn{k2ri)y 
and  the  held  in  region  two  is  given  by 


U^‘^\r,9)  =  By‘'\jyk2r)  cos(n6'). 

n=0 


(4.9b) 


(4.9c) 
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The  superscript  (he)  on  the  expansion  coefficients  in  (4.9)  is  meant  to  remind  the  reader 
that  these  expansion  coefficients  are  for  the  host  cylinder. 

Turning  our  attention  to  the  cylindrical  annulus  shown  in  Figure  4.1(a),  if  we  think  of 
the  cylindrical  annulus  as  the  host  cylinder  into  which  a  scatterer,  namely,  the  cylindrical 
void,  has  been  inserted,  then,  the  total  helds  in  regions  one,  two,  and  three  may  be  written 
as  1/W(r,0)  =  f/W(r,0)  +  =  U^^\r,e)  +  W^^\r,9)  and  V^^\r,e)  = 

9),  respectively.  Noting  that  the  IF-helds  also  satisfy  the  2D  Helmholtz  equation  and 
imposing  the  Sommerfeld  radiation  condition  on  W^^^^\r,9),  as  well  as  requiring  W^^\r,9) 
to  be  regular  at  r  =  0,  yields 

'W^^‘^^\r,9)'\  oo  r  Hnihr) 

W^^\r,9)  J^{k2r)  +  Ct^Yr,{k2r)  cos{n9).  (4.10) 

_W^^\r,9)  \  n=o  ^  Dt^Jnikir) 

The  superscript  (cv)  on  the  expansion  coefficients  in  (4.10)  reminds  us  of  the  presence  of 
the  cylindrical  uoid.  To  hnd  the  unknown  coefficients  in  (4.10),  we  require  that  the  H-£elds 
and  their  normal  derivatives  be  continuous  across  the  r  =  ri  interface,  as  well  as  the  r  =  r2 
interface  to  obtain 


-Hnikiri)  Jn{k2ri)  Yn{k2ri)  0  1 

-H'^{kiri)  KJ'^{k2ri)  KY^{k2ri)  0 

0  Jn{k2r2)  Yn{k2r2)  -Jn{kir2) 

0  KJ^{k2r2)  i^Y^{k2r2)  -Jn{kir2) 

- „ - -  L 

=M 

JnikiTi)  +  Hnikiri)  -  B^n’'^  Jn{k2ri) 

Uk,r,)  +  At^H'^k.r,)  -  kB^^  J'^{k2r,) 

-Blt^^Uk2r2)  ’  ^  ^ 

J'^{k2r2) 

and  the  prime  denotes  the  derivative  with  respect  to  the  argument.  Solving  (4.11)  for  An^^ 
and  using  (4.2b)  to  simplify  the  result,  yields 

det{M)Al^''^  =  -  {Jn+i{kir2)-,  Jn{k2r2)} 

X  {Hn+i{k,ny,Yr,{k2n)}  +  {Jn+i{kiny,Yn{k2n)}) ,  (4.12a) 
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where 


det(M)  =  {Hn+i{kiri)]  Yn{k2ri)}  {Jn+i{kir2)]  Jn{k2r2)} 

-  {Hn+l{kiri)-,  Jn{k2ri)}  {Jn+l{kir2)-,  h^n(/^2»^2)}  • 

To  simplify  (4.12a)  further,  we  use  (4.9b)  and  note  that 


(4.12b) 


{Yn+i{kiri)]  Yn{k2ri)}  {Jn+i{kiri)]  Jn{k2ri)} 

-  {Yn+likiTi)-,  Jn{k2ri)}  {Jn+i{kiri)-,  Yn{k2ri)}  =  ^ 

to  obtain 

Tcv)_  i  (  2  y  {Jn+i(fcir2);Jn(fc2r2)}  + 

dei{M)\7;k,Tj  {H^^,{k,n)-Uk2n)y  ^  ^  -  ^ 

The  coefficients  in  (4.13)  are  the  exact  expansion  coefficients  of  {r,  6) .  To  obtain 

the  approximate  coefficients,  we  hrst  note  that  Yn{k2r2)  ~  — ihfn(^2^"2))  if  |^2k2  ^  1  [13], 

which  allows  us  to  rewrite  (4.13)  as 


^(cv)  ^  i  (JL\  "  1 

”  V^/ciri/  {hf„+i(A;iri);  J„(/C2ri)} 

X 


(4.14) 


{Hn+i{kiri)]Yn{k2ri)}  5n  -  i{Hn+i{kiri)]  Jn{k2ri)}  J  ’ 
where  5n  is  given  by  (4.5).  To  develop  (4.14)  fnrther,  we  note  that  |(5„|  1  if  kir2  1  and 

\k2\r2  ^  1,  as  can  be  seen  from  the  small  argnment  forms  of  Jn{C)  and  Hn{C)  [13]-  Therefore, 
we  can  expand  (4.14)  in  powers  of  (5„  to  hnally  obtain 

2i 


5n, 


n 


e  z+. 


(4.15) 


m 


^Tikin  {Hn+i{kiri)]  Jn{k2ri)} 

Notice  that  the  above  An'''^  coefficients  are  identical  to  the  expansion  coefficients  given 
(4.8)  of  Section  4.3. 

To  numerically  illustrate  the  SCV  approximation,  the  far-held  pattern  of  9)  in 

the  forward  direction,  0  =  0,  as  a  function  of  k2r2  is  show  in  Figure  4.2.  The  far-held  pattern 
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of  9)  is  defined  by 


COSM).  (4.16) 

n=0 

From  Figure  4.2,  we  see  that  the  exact  (computed  with  (4.13))  and  the  approximate  (com¬ 
puted  with  (4.15))  far-held  patterns  are  in  good  agreement  for  small  k2r2,  say,  k2r2  <  0.3. 
Also  from  Figure  4.2,  we  see  that  the  SCV  approximation  becomes  progressively  worse  as 
k2r2  approaches  unity.  This  is  expected  as  the  SCV  approximation  requires  that  both  kir2 
and  k2r2  are  much  smaller  than  unity. 


Figure  4.2:  The  magnitude  and  phase  of  the  far-held  pattern  in  the  forward  direction  for  a 
Tehon  cylindrical  annulus  in  vacuum,  with  an  outer  radius  of  10  cm  at  100  GHz,  is  shown. 
The  permittivity  of  Tehon  at  100 GHz  is  2.05  with  a  negligible  loss-tangent  [10].  In  the 
computation  of  (4.16),  we  only  summed  the  hrst  N  =  \k1r2  +  4(A;ir2)^/^  -|-  2"|  terms  [3, 
Appendix  Cj. 


4.5  Energy  conservation 

In  this  section,  we  present  a  relationship  between  the  rate  at  which  the  energy  is  extin¬ 
guished  by  the  cylindrical  void  from  the  9)  held.  Also,  a  numerical  example  illustrat¬ 

ing  that  the  SCV  approximation  is  in  good  agreement  with  the  derived  energy  conservation 
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relationship  is  presented. 


We  begin  by  constructing  an  imaginary  concentric  cylinder  of  radius  R  >  ri  and  length  L 
around  the  host  cylinder  shown  in  Figure  4.1(b).  Then,  the  rate  at  which  the  energy 

is  absorbed  within  the  imaginary  concentric  cylinder  is  given  by 

/TT 

St/-fd0,  (4.17a) 

•TT 

where  f  =  cos6*x  +  sin^y  (see  Figure  4.1),  and  the  time-averaged  Poynting  vector  is  given 
by 

St' =  (£/'■>)•].  (4.17b) 

In  (4.17b),  Re  denotes  the  real  part,  *  denotes  the  complex  conjugate,  and  k  =  u/c.  Now, 
we  consider  the  rate  at  which  the  energy  is  absorbed  by  the  cylindrical  annulus.  By 

proceeding  as  before,  we  immediately  obtain 


T(abs)  _ 

V  - 


=  -RL  [  Sy  ■  f  d0,  where  Sy  =  -  Re  .  (4.18) 

2  Ank  ^  ^ 


Substituting  R^^^  =  into  (4.18)  and  using  (4.17)  to  simplify  the  result  yields 


^ext  ^  YyFbs)  _  YyFbs) 


(4.19a) 


where 


iR(^)V  ■  f  d0,  (4.19b) 


^(sca)  _ 


i^(sca)y  ('pp(sca)^R  _  ^ 


(4.19c) 


We  interpret 


as  the  rate  at  which  the  energy  is  extinguished  by  a  scatterer,  namely. 


the  cylindrical  void,  in  the  presence  of  the  host  cylinder.  In  other  words,  it  is  the  rate  at 
which  the  energy  is  depleted  by  the  cylindrical  void  from  the  total  held,  outside  the 
host  cylinder.  Moreover,  from  (4.19)  we  see  that  if  the  cylindrical  annulus  is  nonabsorbing 
(e2  is  purely  real)  then  Finally,  substituting  U'^^\r,6)  and  IR*^®'^®')(r,  6*)  (see 
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(4.9a)  and  (4.10))  into  (4.19b),  and  then  integrating  the  result  over  9,  yields 

j-  OO 

w"  =  -  ^  Es"  (R**  K"*]  +  K"'  (4'”’)*] )  ■  (4.20) 

n=0 

We  interpret  the  hrst  term  in  (4.20)  as  the  rate  at  which  extinguishes  energy  from 

and  the  second  term  as  the  rate  at  which  extinguishes  energy  from 

To  illustrate  that  the  SCV  approximation  is  in  good  agreement  with  the  energy  conser¬ 
vation  principle,  we  compute  using  the  exact  and  approximate  coefficients.  Recall 
that  the  exact  coefficients  are  given  by  (4.13),  and  the  approximate  coefficients  by 

(4.15).  The  results  of  the  above-mentioned  computations  are  shown  in  Figure  4.3.  From 
Figure  4.3,  we  see  that  the  SCV  approximation  conserves  energy  to  roughly  10  percent  for 
k2r2  <  1,  which  does  indicate  that  the  SCV  approximation  is  in  good  agreement  with  the 
energy  conservation  principle. 


Figure  4.3:  The  rate  (normalized  by  Lc/Stt)  at  which  energy  is  extinguished  by  the 

cylindrical  void  from  the  total  held  outside  the  host  cylinder  is  shown  as  a  function  of  k2r2- 
The  above  plot  was  produced  with  the  same  parameters  as  the  ones  described  in  the  caption 
of  Figure  4.2. 
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4.6  Conclusions 


In  this  paper,  we  investigated  a  monochromatic  plane  wave  scattering  from  a  solid  ho¬ 
mogeneous  cylinder  (host  cylinder)  and  a  cylindrical  annulus,  referring  to  the  inner  part  of 
the  cylindrical  annulus  as  the  cylindrical  void.  It  was  shown  that  if  the  cylindrical  void  is 
thought  of  as  a  scatterer  inserted  into  the  host  cylinder,  then  the  scattered  field  due  to  the 
cylindrical  void  may  be  approximated  by  the  screened  cylindrical  void  (SCV)  approxima¬ 
tion,  see  Section  4.3.  The  SCV  approximation  was  derived  intuitively  in  Section  4.3  and 
rigorously  in  Section  4.4.  Furthermore,  a  formula  for  the  rate  at  which  energy  is  depleted  by 
the  cylindrical  void  from  the  total  field  outside  the  host  cylinder  was  derived  in  Section  4.5. 
The  numerical  examples  in  Sections  4.4  and  4.5  showed  that  the  SCV  approximation  is  in 
good  agreement  with  the  exact  solution  if  the  cylindrical  void  is  small. 
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CHAPTER  5 


SCATTERING  FROM  A  LARGE  CYLINDER  WITH  AN  ECCENTRICALLY 
EMBEDDED  CORE:  AN  ORDERS-OF-SCATTERING  APPROXIMATION 
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5.1  Abstract 

We  develop  an  orders-of-scattering  approximation,  termed  the  “screened  cylindrical 
void/core”  (SGV)  approximation,  for  a  composite  cylinder.  The  composite  cylinder  consists 
of  a  large  host  cylinder  that  contains  a  small,  eccentrically  embedded,  core  cylinder.  The  SGV 
approximation  is  developed  via  separation  of  variables  in  conjunction  with  addition  theorems 
for  cylindrical  functions.  We  show  that  the  SGV  approximation  is  in  good  agreement  with 
the  numerically-exact  solution.  A  simple  physical  interpretation  of  the  SGV  approximation 
is  also  presented. 

5.2  Introduction 

Gonsider  a  monochromatic  plane  wave  scattering  from  an  infinitely  long  homogeneous 
and  isotropic  composite  cylinder.  The  composite  cylinder  is  composed  of  a  small  core  cylinder 
of  radius  b  that  is  eccentrically  embedded  into  a  large  host  cylinder  of  radius  a,  as  shown  in 
Figure  5.1.  To  experimentally  isolate  the  core  cylinder’s  contribution  to  the  scattered  held  of 
the  composite  cylinder,  one  would  measure  the  total  held  outside  the  composite  cylinder 
and  the  total  held  outside  an  identical  host  cylinder.  Then,  the  difference,  6)  = 
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Figure  5.1:  The  cross-sectional  view  of  the  composite  cylinder,  with  regions  labeled  by  a 
number,  is  shown.  Region  1  is  the  space  outside  of  the  composite  cylinder  (r  >  a).  Region  2 
is  the  host  cylinder,  and  Region  3  is  the  core  cylinder.  The  origin  of  the  (r,  6)  coordinate 
system,  where  —n  <  9  <  n,  is  centered  on  the  host  cylinder,  and  the  origin  of  the  (p,  0) 
coordinate  system,  where  — tt  <  0  <  tt,  is  centered  on  the  core  cylinder.  The  axes  of  these 
two  coordinate  systems  are  parallel  to  each  other  and  the  center  of  the  (p,  0)  coordinate 
system  is  offset  by  vq  cos  x  -|-  tq  sin  9o  y  with  respect  to  the  origin  of  the  (r,  9)  coordinate 
system. 


U^^\r,9)  —  6*),  would  contain  the  effect  that  the  core  cylinder  had  on  the  scattered 

held.  In  our  recent  paper  [1],  we  considered  the  simplest  composite  cylinder  geometry  (the 
core  cylinder  is  concentric  with  the  host  cylinder)  and  developed  an  approximation  to 
which  we  termed  the  “screened  cylindrical  void/core”  (SCV)  approximation.  In  this  paper, 
we  derive  an  analogous  formula  for  an  eccentrically  stratihed  composite  cylinder,  which  can 
also  be  interpreted  as  an  orders-of-scattering  approximation.  Furthermore,  we  numerically 
investigate  the  accnracy  of  the  SCV  approximation  when  \k2\a  ~  300  and  0  <  \k3\b  <  1, 
where  k2  {k^)  is  the  wavennmber  in  the  host  (core)  cylinder. 

Scattering  by  an  eccentrically  stratihed  composite  cylinder  has  previously  been  considered 
in  the  literature  in  various  contexts  [2-4]  and  by  various  techniques  [5-7] .  In  the  electromag¬ 
netic  context,  a  pertnrbation  series  solution  has  been  constructed  in  powers  of  (/C3  — /C2)  [8,  9], 
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b  [10],  and  eccentricity  [11,  12]  by  using  separation  of  variables.  An  “exact”  treatment  based 
on  separation  of  variables  with  a  truncation  of  the  resultant  inhnite  size  matrix  is  also  avail¬ 
able,  e.g.  in  [13].  Our  orders-of-scattering  approach  is  also  based  on  separation  of  variables, 
but  the  resultant  power  series  expansion  of  the  solution  is  different  from  the  ones  mentioned 
above. 

There  are  many  diverse  applications  where  the  scattering  by  an  eccentrically  stratihed 
composite  cylinder  is  important;  for  example,  see  [4,  8,  9]  and  references  therein.  As  men¬ 
tioned  in  [1] ,  we  are  particularly  interested  in  using  the  composite  cylinder  to  study  Anderson 
localization  [14-16]  at  millimeter/sub- millimeter  wavelengths.  At  these  wavelengths,  both 
the  amplitude  and  the  phase  of  the  electromagnetic  held  can  be  easily  measured  with  a 
vector  network  analyzer  [17],  and  the  preparation  of  disordered  samples  is  straightforward 
with  standard  computer-numerically-controlled  milling  techniques.  For  example,  a  sample 
may  be  prepared  by  drilling  small  holes  in  a  large  Tehon  (ultra  low-loss  material)  cylinder. 
Furthermore,  essentially  arbitrary  realizations  of  the  same  random  disorder  may  be  gener¬ 
ated  by  putting  the  sample  on  a  rotational  stage  and  illuminating  it  from  the  side.  When 
the  number  of  small  scatterers  is  large,  say,  over  a  thousand,  then  what  is  important  is  the 
rate  at  which  the  scatterer  extinguishes  the  energy  from  the  incident  held,  rather  than  the 
geometrical  shape/size  of  each  individual  scatterer  [18,  19].  Practically,  the  host  cylinder 
needs  to  be  rather  large  (a  ~  10  cm)  in  order  to  accommodate  thousands  of  small  holes 
{b  ~  0.3  mm);  hence  the  numerical  examples  considered  in  this  paper  are  for  k2a  ~  300  and 
0  <  k^b  <  1.  Therefore,  the  physical  insight  gained  from  considering  the  scattering  by  a 
single  core  cylinder  eccentrically  embedded  into  a  large  host,  as  discussed  in  this  paper,  may 
be  of  beneht  in  understanding  the  experimental  model  described  above.  Also,  the  approach 
taken  in  this  paper  may  pave  the  way  for  methods  that  may  accurately  describe  the  full 
envisioned  experiment,  where  thousands  of  core  cylinders  are  eccentrically  embedded  into 
one  large  host  cylinder. 
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5.3  Background  and  conventions 


Throughout  the  paper,  we  will  use  the  Gaussian  unit  system,  and  we  will  assume  that  all 
helds  are  harmonic  in  time  with  a  exp(— lent)  time  factor,  where  u  is  the  angular  frequency. 
Furthermore,  we  will  assume  that  all  helds  are  polarized  in  the  positive  z-direction.  The 
positive  z-direction  is  out  of  the  page  in  Figure  5.1.  All  media  considered  in  this  paper 
are  assumed  to  be  non-magnetic,  and  the  permittivity  of  the  host  and  the  core  cylinder  are 
denoted  by  e-2  and  63,  respectively.  Furthermore,  ei  denotes  the  permittivity  of  the  space 
outside  the  composite  cylinder  (Region  1  in  Figure  5.1)  and  is  assumed  to  be  purely  real, 
whereas  62  and  63  may  be  complex. 

Let  us  note  that  all  helds  in  this  paper  satisfy  the  two-dimensional  (2D)  Helmholtz  equa¬ 
tion.  The  radial  solution  of  the  2D  Helmholtz  equation  is  composed  of  a  linear  combination 
of  an  integer  order  Bessel  function  of  the  hrst  kind  and  an  integer  order  Hankel  function  of 
the  hrst  kind,  which  we  denote  by  Jm{,i)  and  respectively.  The  functions  Jm(0  and 

satisfy  the  Wronskian  relationship  [20,  §9.1] 


UOKSi)  -  J'jonuo 


and  the  recurrence  relation  [20,  §9.2] 


2i 


(5.1a) 
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«;„(«)  =  -  'tm+iK).  (5-lb) 

where  T  denotes  J  or  H,  and  the  prime  denotes  the  derivative  with  respect  to  the  argument. 
It  is  convenient  to  introduce  the  shorthand  curly  bracket  notation,  {Tm+i(0j ‘^m(h)})  by 
which  we  mean 


{'Lm+l(0)  *^m(h)}  —  ^m+1  ('C)*^m(h)  ^  ^m(0*^m+l  (h)- 

For  example,  if  T  and  <h  satisfy  (5.1b),  then 

(5.1c) 


80 


Lastly,  we  note  the  Jacobi-Anger  expansion  of  a  plane  wave  [21,  p.  37],  namely, 

ei«c°«®  =  (5.2) 

m 

where  ^  indicates  the  summation  from  m  =  —oo  to  m  =  oo. 

m 

5.4  Host  cylinder 


Consider  a  unit  plane  wave,  0^'“^  =  exp  {ikir  cos  0),  incident  from  Region  1  onto  the  host 
cylinder,  see  Figure  5.1  with  b  =  0  (i.e.,  without  the  core  cylinder).  Then,  after  decomposing 
the  total  held  in  Region  1  as  we  have  [1] 


■Uhca)(^,^)- 

U(2)(r,0) 


AmO 


(5.3) 


where  ki  =  ^iUjjc  for  i  =  1,  2  and  c  is  the  speed  of  light  in  vacuum.  In  (5.3),  denotes 
the  total  held  inside  the  host  cylinder,  and  the  expansion  coefficients  are  given  by 


{'7m+l(^l®))  Jmik20j')} 
{.^m+1  (^1®) )  '7m(^2®)} 


'Kkia{Hm+i{,kia)-,Jm{,k2a)}' 


(5.4a) 

(5.4b) 


5.5  Composite  cylinder 


If  the  plane  wave  is  incident  from  Region  1  onto  the  composite  cylinder  shown  in 

Figure  5.1,  then  the  total  helds  in  Regions  1,  2,  and  3  may  be  written  as 


and 


respectively,  where 


[/(i)  (r,  9)  =  (r,  9)  +  (r,  9) , 

(5.5a) 

7/(2)  (r,  9- p,ct))  =  U(2)  (r,  9)  +  0)  ^ 

(5.5b) 

m 

(5.5c) 
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(5.5d) 


m 

e-  P,  0)  =  i™  +  C^Hyk^py^^)  ,  (5.5e) 

m 

and  ks  =  yBibjjc.  In  writing  (5.5),  we  are  thinking  of  the  composite  cylinder  as  the  host 
cylinder  into  which  a  scatterer  (the  core  cylinder)  has  been  inserted.  Also,  notice  that  we 
required  U^^\p^(f))  to  be  finite  at  p  =  0,  and  imposed  the  Sommerfeld  radiation  {outgoing 
cylindrical  wave)  condition  on  V^‘^^^^{r,9).  To  hnd  the  unknown  expansion  coefficients  in 
(5.5),  we  require  that  the  electric  field  and  its  normal  derivative  be  continuous  across  the 
p  =  b  and  r  =  a  interfaces. 

To  apply  the  continuity  conditions  at  the  p  =  b  interface,  we  hrst  re-express  9;  p,  0) 

solely  in  terms  of  the  (p,  0)  coordinate  system  by  using  Graf’s  addition  theorem  [21,  §2.5], 
[20,  §9.2].  Namely,  using 

Jm{k2ry^^  =  J„^-n(A:2ro)e'(™-")^°J4A;2p)e’"^ 

n 

and  (5.3)  with  (5.2),  we  obtain 

y^\p,  y  =  i"”^(^2P)e'"^  (B^  +  By +  J2  rCmHyk2py^^, 

n  m  m 

where  Tnm  =  Then,  requiring  that  and 

on  p  =  6  yields 

DpJp{ky)  =  CpHp{k2b)  +  {-lYJp{k2b)  ^  Tpm  (B,„  -h  By  ,  (5.6a) 

m 

and 

^DpJ'YkY)  =  CpHYk2b)  +  {-lYJ'p{k2b)y^Tpm  (B™  +  By  ,  (5.6b) 

respectively.  Eliminating  Dp  from  (5.6)  yields 

Cp  =  (— l)^Ap  Tpm  (Bm  -|-  Bm) ,  (5.7a) 
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where 


A  _  {»^p+l(^3^);  Jp{k2h)"\ 

{jp^^{hby,Hp{hb)y  ’ 

Similarly,  to  apply  the  continuity  conditions  at  the  r  =  a  interface,  we  first  re-express 
(r,  6]  p,  0)  solely  in  terms  of  the  (r,  9)  coordinate  system  by  using  Graf’s  addition  theorem 
for  Hm{k2p)e^"^‘^  [21,  §2.5],  [20,  §9.2],  Namely,  using 

n 

for  r  >  ro,  and  (5.3)  with  (5.2),  we  obtain 

m  n  m 

ioi  tq  +  b  <  r  <  a.  Then,  requiring  that  and  on  r  =  6  yields 

(Bp  -|-  Bp)  Jp{k2a)  -\-  Hp{k2Ci)  ^^(— l)™'TpmGm  =  Jp{kia)  +  (Ap  -|-  Ap)  Hp(kia),  (5.8a) 

m 

and 

(Bp  +  Bp)  J'p{k2a)  +  H'p{k2a)  V(-l)”^Tp^G™  =  ^  [j;(A;ia)  +  (Ap  +  Ap)  H'p{kra)]  ,  (5.8b) 

m 

respectively.  To  solve  (5.8)  for  Ap  in  terms  of  Cm,  we  eliminate  (Bp  -|-  Bp)  from  (5.8),  and 
then  use  (5.1a)  and  (5.4)  to  rewrite  the  result  as 

Ap  =  Bp5^(-l)™Tp,„G^.  (5.9) 

m 

To  solve  (5.8)  for  Bp  in  terms  of  Ap,  we  eliminate  Cm  from  (5.8),  and  use  (5.1a)  to  obtain 


2i 


+  Bp)  —  (Ap  -|-  Ap)  {ifp+i(A;ia);  i/p(/c2a)}  +  {Jp+i{kia)]  Hp{k2a)}  .  (5.10) 


T  \^p 

nkia 

To  simplify  (5.10)  further,  we  substitute  (5.4)  into  (5.10)  and  note  that 
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"I "^p+i (^1®) j  Hp(^k2Ci)}  {Hp+ii^kici^ ^  Jp(^k2Ci)} 

{  Jp+i (^1®) )  '-^p(^2®) }  {Hp+i i^kici^ ,  Hp(^k2Ci) } 


to  obtain 

Bp  (/uio) ,  Hp(^k2(i)}  -Ap. 

Finally,  substitnting  (5.11)  into  (5.7a),  and  pntting  the  resnlt  into  (5.9)  yields 

^  ^  i^^mn  Fmn) 

n 

where 


Fmn  =  ^  T^p ApTp„  {i7„+i (/cio) ;  Hn{k2a)}  , 


Gm  —  Tjnp^pTpn  j  Bn; 

n  \  p  / 


and  drrm  denotes  the  Kronecker  delta  fnnction. 


(5.11) 

(5.12a) 

(5.12b) 

(5.12c) 


Notice  that  in  (5.12)  the  core  cylinder  parameters,  namely,  k^  and  b,  are  solely  contained 
in  Ap,  see  (5.7b).  Fnrthermore,  from  (5.7b)  and  the  small  argument  forms  of  Jp  and  Hp,  we 
see  that  if  the  core  cylinder  is  small,  then  so  is  Ap.  This  suggests  that  (5.12a)  can  be  solved 
via  the  Neumann  series  (Taylor  series  expansion,  if  you  will),  i.e.. 


A  =  (I-F)"^G  =  ^F^G,  (5.13) 

i=0 

where  Fmn,  Gm  are  the  elements  of  A,  F,  G,  respectively,  and  I  is  the  identity  matrix. 
The  Neumann  series  in  (5.13)  converges,  provided  that  the  spectral  radius  of  F  is  less  than 
one  [22,  §4.3].  The  spectral  radius  of  F  for  a  large  host  cylinder,  \k2\a  ~  300,  with  an 
eccentrically  embedded  core  cylinder  is  shown  in  Figure  5.2.  From  Figure  5.2,  we  see  that 


the  spectral  radius  of  F  is  indeed  much  smaller  than  one  and  thus,  we  expect  the  Neumann 
series  in  (5.13)  to  converge  rapidly.  We  will  discuss  the  spectral  radius  of  F  further  in  Sec.  5.7, 
but  for  now  turn  our  attention  to  the  physical  interpretation  of  the  SCV  approximation. 
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Figure  5.2:  The  spectral  radius  of  F  at  100  GHz  for  a  Teflon  host  cylinder  (a  =  10  cm) 
with  an  eccentrically  embedded  quartz  core  cylinder  is  shown  as  a  function  of  and 

eccentricity,  r^/a  (with  =  0).  The  permittivity  of  Teflon  and  quartz  at  100  GHz  is  2.1  and 
3.8  with  a  negligible  loss-tangent  [23],  respectively. 


5.6  The  SCV  approximation  and  its  physical  interpretation 

If  only  the  £  =  0  term  is  retained  in  (5.13),  we  obtain  the  SGV  approximation,  namely, 

—  Gm  =  j  ®n-  (5-14) 

n  \  p  / 

To  interpret  (5.14)  physically,  we  consider  the  following  three-step  scattering  process: 

1.  If  a  unit  plane  wave,  =  exp  (ifcir  cos^),  is  incident  on  the  host  cylinder,  then  the 
field  inside  the  host  cylinder,  I[J*^^^(r,  0),  is  given  by  (5.3).  Rewriting  terms  of  the 
(p,  0)  coordinate  system  yields 

U(2)  (p,0)  =  (5.15) 
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2.  If  we  use  (5.15)  as  an  incident  field  for  the  core  cylinder,  then  the  resulting  scattered 


field  is 


(5.16a) 

m 

and  the  held  inside  the  core  cylinder  is 

U^^\pA)  =  (5.16b) 

m 

Substituting  (5.15)  and  (5.16)  into  the  continuity  conditions  for  the  p  =  b  interface, 
and  eliminating  from  the  resultant  two  equations,  yields 

m 

3.  Finally,  if  we  use  (5.16a)  with  (5.17)  as  an  incident  held  (from  within  the  host  cylinder) 
on  the  r  =  a  interface,  then  there  will  be  an  outgoing  held  outside  the  host  cylinder 
given  by 


y(sca)(^,^)  _  (5.18a) 

m 

and  a  regular  (hnite  at  r  =  0)  held  inside  the  host  cylinder  given  by 


J2^^BmJm{k2r)e^^^.  (5.18b) 

m 

Rewriting  (5.16a)  in  terms  of  the  {r,6)  coordinate  system  and  substituting  it,  as  well 
as  (5.18),  into  the  continuity  conditions  for  the  r  =  a  interface,  and  eliminating  Bm 
from  the  resultant  two  equations  yields 


TmpApTpn  j 

n  \  p  / 


(5.19) 


By  comparing  (5.19)  with  (5.14),  we  conclude  that  the  SCV  approximation  can  be  viewed 
as  an  orders-of-scattering  approximation.  Moreover,  from  the  above  three-step  scattering 
process,  we  see  that  is  the  “screening”  ehect  of  the  host  cylinder  on  and 

is  the  “screening”  ehect  of  the  host  cylinder  on  These  two  screening  ehects  are  identical 


if  the  core  cylinder  is  concentric  with  the  host  cylinder,  as  we  have  shown  in  [1],  To  see 
that  (5.19),  or  equivalently  (5.14),  reduces  to  our  previous  result,  we  note  that  Tpn  = 
and  Tmp  =  when  ro  =  0,  and  thus,  the  sums  in  (5.19)  collapse  and  we  obtain 

Am  =  B^A^. 

5.7  Numerical  examples  and  limitations 

In  practice,  the  computation  of  the  Am  coefficients  via  (5.12)  or  (5.14)  requires  the 
truncation  of  the  inhnite  sums,  as  well  as  the  index  m.  From  (5.12b),  (5.12c),  and  (5.14), 
we  see  that  the  sum  over  p  is  controlled  by  the  small  core  cylinder  parameters,  namely,  Ap. 
This  observation  suggests  that  the  summation  over  p  be  terminated  at  pmax  (i-e.,  |p|  <  Pmax), 
where  Pmax  is  given  by  the  well-known  Wiscombe’s  criterion  for  small  scatterers  [24],  namely. 


Pmax 


k2b  +  4{k2bf^  +  1 


(5.20a) 


The  sum  over  n,  as  well  as  the  index  m,  are  controlled  by  the  large  host  cylinder  and  thus, 
they  are  terminated  at  iVmax  (i.e.,  \n\  <  N^ax  and  \m\  <  N^ax),  where  N^ax  is  given  by  the 
Wiscombe’s  criterion  for  relatively  large  scatterers  [24],  namely. 


N  = 

’  mfl.x 


kia  +  4.05  {kio)^^^  +  2 


(5.20b) 


We  note  that  a  termination  criterion  in  terms  of  prescribed  relative  error  has  become  available 
recently  [25],  but  for  our  purposes,  the  termination  condition  given  by  (5.20)  will  be  sufficient. 

To  numerically  illustrate  the  accuracy  of  the  SCV  approximation,  we  compute  the  relative 
error  in  the  rate  at  which  the  energy  is  extinguished  by  the  core  cylinder  in  the  presence  of 
the  host  cylinder.  The  rate  at  which  the  energy  (per  unit  length  of  the  composite  cylinder) 
is  depleted  by  the  core  cylinder  from  the  total  field,  outside  the  the  host  cylinder  is 
given  by  [1] 

2  -^max 

(Re[A,„]  +  2Re[A,„A^*]),  (5.21) 

717=  A^max 
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where  Re  denotes  the  real  part  and  *  denotes  the  complex  conjugate.  We  compute  the 
SCV  approximate  and  numerically-exact  (~  7  significant  digits)  by  using  (5.21)  with 
(5.14)  and  (5.21)  with  (5.13),  respectively.  The  top  row  of  Figure  5.3  shows  that  the  SCV 
approximation  is  in  good  agreement  with  the  numerically-exact  solution,  and  the  bottom 
row  of  Figure  5.3  demonstrates  that  the  Neumann  series  in  (5.13)  converges  rapidly  as  one 
would  expect  from  the  spectral  radius  of  F,  see  Figure  5.2.  Furthermore,  from  Figure  5.3  we 
see  that  the  relative  error  in  is  almost  independent  of  the  angular  position  of  the  core 
cylinder  but  does  depend  on  its  radial  position,  see  Figure  5.3  with  ro/a  >  0.7. 

The  dependence  of  the  relative  error  in  on  the  radial  position  of  the  core  cylinder 
may  be  explained  in  terms  of  the  internal  resonances  of  the  host  cylinder.  These  resonances 
are  often  referred  to  as  Mie  resonances,  morphological  resonances,  whispering-gallery  modes, 
or  natural/eigen  modes.  At  100  GHz,  the  10  cm  host  cylinder  is  about  hundred  times  larger 
than  the  wavelength  of  the  incident  light  and  thus,  the  interaction  of  light  with  the  host 
cylinder  can  be  described  by  ray  theory.  If  a  ray  inside  the  host  cylinder  strikes  the  surface 
of  the  host  cylinder  above  the  critical  angle,  then  the  ray’s  trajectory  will  be  bounded  by  a 
cylindrical  annulus  with  outer  radius  a  and  inner  radius  rcaustic-  To  find  the  caustic  radius, 
Raustic,  we  set  the  ray’s  angular  momentum  \k2fi\rh  equal  to  \m\h  (the  angular  momentum 
of  the  mth  eigenmode)  and  note  that  /c|  =  ^  +  /c2  r  obtain 


^caustic 


m 

k2 


(5.22a) 


In  the  derivation  of  (5.22a),  we  used  the  fact  that  the  radial  component  of  the  wavevector 
must  vanish  on  rcaustic,  he.,  k2,r{j'  =  ''"caustic)  =  0  [26,  27].  Furthermore,  we  can  deduce  the 
range  of  potentially  excited  eigenmodes  of  the  host  cylinder  as  follows.  If  a  ray  inside  the  host 
cylinder  strikes  the  surface  at  an  angle  7  with  respect  to  the  normal,  then  by  equating  the 
ray’s  and  modal  angular  momenta  (|m|  =  |A;2|asin7),  and  using  the  total  internal  reflection 


condition,  a/ 61/62  <  sin  7  <  1,  we  obtain 

|A;i|a  <m<  \k2\a. 


(5.22b) 
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Figure  5.3:  The  relative  error  in  (in  percent)  is  shown  as  a  function  of  \k3\b  and  eccen¬ 
tricity,  ro/a,  for  various  60  angles.  The  top  row  shows  the  relative  error  if  only  the  £  =  0  term 
is  retained  in  (5.13),  i.e.,  the  SCV  approximation,  and  the  bottom  row  shows  the  relative 
error  if  the  £  =  0  and  £  =  1  terms  are  retained.  The  above  plot  was  produced  with  the  same 
parameters  as  the  ones  described  in  the  caption  of  Figure  5.2 


Finally,  from  (5.22),  we  see  why  the  SCV  approximation  worsens  when  the  radial  location 
of  the  core  cylinder  exceeds  the  caustic  radius,  see  Figure  5.3  for  ro/a  >  Ucaustic/a  ~  0.7. 

If  the  frequency  of  the  incident  wave  corresponds  to  one  of  the  eigenfrequencies  of  the  host 
cylinder,  then  the  Neumann  series  in  (5.13)  will  fail  to  converge  only  when  tq  >  Ucaustic-  For 
example,  the  mode  m  =  228  is  excited  in  resonance  at  approximately  99.823859  GHz,  i.e., 
the  denominator  of  B228  vanishes  at  this  frequency and  the  spectral  radius  of  F  exceeds 
unity  when  ro/a  <  rcaustic/a  =  228/ {k2a)  ~  0.75  as  shown  in  Figure  5.4.  Moreover,  from 
Figure  5.2  we  see  that  the  SCV  approximation  remains  valid  even  at  resonance  frequency, 
provided  that  ro/a  <  rcaustic/a  ~  0.75. 


Figure  5.4:  The  spectral  radius  of  F  at  eigenfrequency  99.823859  GHz  is  shown  as  a  function 
1^3 1 6  and  eccentricity,  ro/a  (with  9o  =  0).  The  above  plot  was  produced  with  the  same 
parameters  as  the  ones  described  in  the  caption  of  Figure  5.2 


^^Strictly  speaking,  this  occurs  at  a  complex  eigenfrequency,  where  the  imaginary  part  of  the  eigenfrequency 
is  related  to  the  spectral  width  of  the  mode  [28]. 


90 


5.8  Conclusions 


In  this  paper,  we  have  extended  the  screen  cylindrical  void/core  (SCV)  approximation  [1] 
to  a  case  where  the  small  core  cylinder  is  eccentrically  embedded  into  a  large  host  cylinder. 
We  physically  interpreted  the  SCV  approximation  as  the  screening  effect  of  the  host  cylinder 
on  the  incident  plane  wave  and  the  wave  scattered  by  the  core  cylinder  (see  Section  5.6). 
Furthermore,  we  showed  that  the  SCV  approximation  may  be  thought  of  as  an  orders-of- 
scattering  approximation. 

The  accuracy  of  the  SCV  approximation  was  demonstrated  numerically  for  an  envisioned 
localization  experiment,  where  a  large  host  cylinder  {k2a  ~  300)  contains  a  small  {k^b  ~  1) 
eccentrically  embedded  core  cylinder.  In  general,  the  SCV  approximation  was  shown  to  be 
in  good  agreement  with  the  exact  solution,  even  at  the  eigenfrequencies  of  the  host  cylinder. 
We  showed  that  if  the  incident  frequency  corresponds  to  one  of  the  eigenfrequencies  of  the 
host  cylinder,  then  the  SCV  approximation  remains  valid,  provided  that  the  eccentricity 
ro/a  does  not  exceed  the  caustic  radius  of  the  mode  (see  Section  5.7).  This  condition  was 
derived  by  considering  the  interplay  between  the  ray  and  wave  pictures  of  the  scattering 
process.  Moreover,  the  ray  picture  offered  a  valuable  physical  insight  into  the  validity  of  the 
SCV  approximation. 

5.9  Acknowledgment 

This  material  is  based  upon  work  supported  in  part  by  the  U.S.  Office  of  Naval  Research  as 
a  Multi-disciplinary  University  Research  Initiative  on  Sound  and  Electromagnetic  Interacting 
Waves  under  grant  number  N00014-10-1-0958. 

5.10  References  Cited 

[1]  A.  J.  Yuffa,  J.  A.  Scales,  Measuring  the  void:  Theoretical  study  of  scat¬ 
tering  by  a  cylindrical  annulus,  J.  Quant.  Spectrosc.  Radiat.  Transfer  (2013). 
doi;10.1016/j.jqsrt.2013.02.017. 


91 


[2]  J.  A.  Roumeliotis,  N.  B.  Kakogiannos,  Acoustic  scattering  from  an  infinite  cylinder  of 
small  radius  coated  by  a  penetrable  one,  Journal  of  the  Acoustical  Society  of  America 
97  (1995)  2074-2081. 

[3]  J.-M.  A.  Noel,  A.  J.  Patitsas,  Modes  of  vibration  in  an  ideal  fluid  bounded  by  two 
eccentric  rigid  infinite  circular  cylindrical  boundaries,  Canadian  Journal  of  Physics  76 
(1998)  729-738. 

[4]  L.-W.  Cai,  Scattering  of  elastic  anti-plane  shear  waves  by  mnltilayered  eccentric  scat- 
terers,  Qnarterly  Jonrnal  of  Mechanics  and  Applied  Mathematics  58  (2005)  165-183. 

[5]  R.  P.  Shaw,  T.  George,  Time  harmonic  aconstic  radiation  from  nonconcentric  circnlar 
cylinders.  The  Jonrnal  of  the  Acoustical  Society  of  America  56  (1974)  1437-1443. 

[6]  X.  Yuan,  D.  R.  Lynch,  J.  W.  Strohbehn,  Conpling  of  hnite  element  and  moment  methods 
for  electromagnetic  scattering  from  inhomogeneous  objects.  Antennas  and  Propagation, 
IEEE  Transactions  on  38  (1990)  386-393. 

[7]  E.  B.  Danila,  J.  M.  Conoir,  J.  L.  Izbicki,  Generalized  Debye  series  expansion:  Part 
II,  treatment  of  eccentric  flnid-solid  cylindrical  interfaces,  Acta  Acnstica  nnited  with 
Acustica  84  (1998)  38-44. 

[8]  N.  K.  Uznnogln,  J.  G.  Fikioris,  Scattering  from  an  inhnite  dielectric  cylinder  embedded 
into  another.  Journal  of  Physics  A:  Mathematical  and  General  12  (1979)  825-834. 

[9]  G.  A.  Valagiannopoulos,  Electromagnetic  scattering  from  two  eccentric  metamaterial 
cylinders  with  freqnency- dependent  permittivities  differing  slightly  each  other.  Progress 
In  Electromagnetics  Research  B  3  (2008)  23-34. 

[10]  J.  A.  Ronmeliotis,  N.  B.  Kakogiannos,  Scattering  from  an  inhnite  cylinder  of  small 
radius  embedded  into  a  dielectric  one.  Microwave  Theory  and  Techniqnes,  IEEE  Trans¬ 
actions  on  42  (1994)  463-470. 

[11]  J.  A.  Ronmeliotis,  J.  G.  Fikioris,  G.  P.  Gonnaris,  Electromagnetic  scattering  from  an 
eccentrically  coated  inhnite  metallic  cylinder,  Jonrnal  of  Applied  Physics  51  (1980) 
4488-4493. 

[12]  G.  P.  Zonros,  J.  A.  Ronmeliotis,  G.-T.  Stathis,  Electromagnetic  scattering  by  an  inhnite 
cylinder  of  material  or  metamaterial  coating  eccentrically  a  dielectric  cylinder,  J.  Opt. 
Soc.  Am.  A  28  (2011)  1076-1085. 

[13]  R.  P.  Parrikar,  A.  A.  Kishk,  A.  Z.  Elsherbeni,  Scattering  from  an  impedance  cylinder 
embedded  in  a  nonconcentric  dielectric  cylinder.  Microwaves,  Antennas  and  Propaga¬ 
tion,  lEE  Proceedings  H  138  (1991)  169-175. 


92 


[14]  P.  W.  Anderson,  Absence  of  diffusion  in  certain  random  lattices,  Phys.  Rev.  109  (1958) 
1492-1505. 

[15]  A.  Laffendiik,  B.  van  Tiggelen,  D.  S.  Wiersma,  Fifty  years  of  Anderson  localization. 
Physics  Today  62  (2009)  24-29. 

[16]  T.  Sperling,  W.  Biihrer,  C.  M.  Aegerter,  G.  Maret,  Direct  determination  of  the  transi¬ 
tion  to  localization  of  light  in  three  dimensions.  Nature  Photonics  7  (2013)  48-52. 

[17]  J.  A.  Scales,  L.  D.  Carr,  D.  B.  McIntosh,  V.  Freilikher,  Y.  P.  Bliokh,  Millimeter  wave 
localization:  Slow  light  and  enhanced  absorption  in  random  dielectric  media,  Phys. 
Rev.  B  76  (2007)  085118. 

[18]  M.  Rusek,  A.  Orlowski,  Analytical  approach  to  localization  of  electromagnetic  waves  in 
two-dimensional  random  media,  Phys.  Rev.  E  51  (1995)  R2763-R2766. 

[19]  M.  Rusek,  A.  Orlowski,  Example  of  self-averaging  in  three  dimensions:  Anderson  local¬ 
ization  of  electromagnetic  waves  in  random  distributions  of  pointlike  scatterers,  Phys. 
Rev.  E  56  (1997)  6090-6094. 

[20]  M.  Abramowitz,  I.  A.  Stegun  (Eds.),  Handbook  of  Mathematical  Functions,  Dover,  New 
York,  1965. 

[21]  P.  A.  Martin,  Multiple  Scattering,  Cambridge  University  Press,  Cambridge,  2006. 

[22]  N.  J.  Higham,  Functions  of  Matrices,  SIAM,  Philadelphia,  2008. 

[23]  P.  F.  Goldsmith,  Quasioptical  Systems:  Gaussian  Beam  Quasioptical  Propagation  and 
Applications,  Wiley-IEEE  Press,  New  York,  1998. 

[24]  W.  J.  Wiscombe,  Improved  Mie  scattering  algorithms,  Appl.  Opt.  19  (1980)  1505-1509. 

[25]  A.  A.  R.  Neves,  D.  Pisignano,  Effect  of  hnite  terms  on  the  truncation  error  of  Mie 
series.  Opt.  Lett.  37  (2012)  2418-2420. 

[26]  G.  Roll,  G.  Schweiger,  Geometrical  optics  model  of  Mie  resonances,  J.  Opt.  Soc.  Am. 
A  17  (2000)  1301-1311. 

[27]  J.  Schulte,  G.  Schweiger,  Resonant  inelastic  scattering  by  use  of  geometrical  optics,  J. 
Opt.  Soc.  Am.  A  20  (2003)  317-324. 

[28]  P.  R.  Conwel,  P.  W.  Barber,  C.  K.  Rushforth,  Resonant  spectra  of  dielectric  spheres, 
J.  Opt.  Soc.  Am.  A  1  (1984)  62-67. 


93 


CHAPTER  6 


GENERALIZATION  OF  THE  SCV  APPROXIMATION  TO  A  CLUSTER  OF 
ECCENTRICALLY  EMBEDDED  CORES 

The  SCV  approximation  derived  for  a  composite  cylinder  with  one  eccentrically  embedded 
core  cylinder  may  be  generalized  to  a  case  of  N  eccentrically  embedded  core  cylinders  via 
the  so-called  cluster  T-matrix.  The  cluster  T-matrix  [1,  2]  [3,  §5.9]  [4,  §7.11]  relates  the 
expansion  coefficients  of  the  wave  incident  on  the  cluster  to  the  expansion  coefficients  of  the 
wave  scattered  by  the  cluster.  In  other  words,  the  cluster  T-matrix  allows  us  to  treat  all 
of  the  eccentrically  embedded  cylinders  as  one  unit.  Therefore,  the  derivation  of  the  SCV 
approximation  for  N  eccentrically  embedded  core  cylinders  will  closely  parallel  the  derivation 
in  Chapter  5. 

The  composite  cylinder  with  N  eccentrically  embedded  core  cylinders  is  shown  in  Fig¬ 
ure  6.1,  where  a  is  the  radius  of  the  host  cylinder  and  the  radius  of  the  ith  core  cylinder  is 
denoted  by  h.  In  what  follows,  it  is  convenient  to  let 


^nikjv)  =  Hn{kjr)e''^\  (6.1a) 

^n{kjr)  =  (6.1b) 

where  kj  is  the  wavenumber  in  the  jth  region,  and  express  Graf’s  addition  theorems  [4,  §2.5] 
[5,  §10.23]  in  the  notation  of  (6.1)  as 


T^(/cw)  =  y^T^_„(/cu)T^(fcv), 

n 

'^raikw)  =  T^_„(-fcu)T^(/cw), 

n 

'^mikw)  =  T^_„(-fcu)T^(fcw), 


(6.2a) 

w  <  u  , 

(6.2b) 

w  >  u  . 

(6.2c) 
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Figure  6.1:  The  cross-sectional  view  of  the  composite  cylinder  is  shown.  Region  iV  -|-  2  is 
the  space  outside  of  the  composite  cylinder  (r  >  a),  region  iV  -|-  1  is  the  host  cylinder,  and 
regions  1, ...  ,7V  are  the  core  cylinders.  The  origin  of  the  global  {r,6)  coordinate  system, 
where  — vr  <  6^  <  tt,  is  centered  on  the  host  cylinder,  and  the  origin  of  the  local  (p*,  0*) 
coordinate  system,  where  — tt  <  4>i  <  tt,  is  centered  on  the  7th  core  cylinder.  The  axes  of  the 
(r,  9)  and  (pj,  coordinate  systems  are  parallel  to  each  other  and  the  center  of  the  (p*,  0j) 
coordinate  system  is  offset  by  r*  with  respect  to  the  origin  of  the  (r,  6)  coordinate  system. 

6.1  Host  cylinder  with  and  without  core  cylinders 

If  we  assume  that  the  incident  wave  on  the  host  cylinder  is  regular  (non-singular)  and 
satisfies  the  2D  Helmholtz  equation,  then  the  total  held  outside  and  inside  the  host  cylinder 
is  given  by 


U(^+2)(r) 

=  y^A^T^(fcjv+2r)  +Bm^{kN+2r), 

(6.3a) 

U(^+i)(r) 

=  y^B^T(fcAf+ir), 

(6.3b) 

m 


respectively.  The  expansion  coefficients  in  (6.3)  are  given  by 
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_  _  {Jm+l{kN+20)  7  Jm  (/cAT+ia)} 
{H^+i{kN+2a)]  Jm{kN+ia)] 

-2i  ^ 

TlkM+2(^{Hm+l{k]\}j^20)]  Jm{kM+lO)} 


(6.4a) 

(6.4b) 


where  the  coefficients  are  assumed  to  be  the  known  expansion  coefficients  of  the  incident 
wave.  If  we  insert  the  cluster  of  core  cylinders  into  the  host  cylinder  in  the  presence  of  the 
U(^+2)  ^(Af+i)  fields,  then  the  new  fields  are  given  by 


[7(^+2)  (r)  =  u(^+2)(r)  +  (6.5a) 

=  I[J*^^+^)(r)  +  I/*^'^’''^^(r),  Tmax  <r<a,  (6.5b) 


where 


=  Y,Am^m{kN+2r), 

m 

f/CA'+b^r)  =  ^  j5^$^(/cjv+ir)  +  C^'^rnikN+ir)  , 


(6.5c) 

(6.5d) 


and  Tmax  is  the  radius  of  the  smallest  imaginary  cylinder  that  circumscribes  the  cluster  of  core 
cylinders,  see  Figure  6.1.  It  is  important  to  stress  that  the  <  r  <  a  condition  in  (6.5b) 
allows  us  to  treat  the  cluster  of  core  cylinders  as  one  entity.  Physically,  we  may  interpret 
U(^+i)(r)  +  E™5A(A:iv+ir)  as  the  incident  held  on  the  cluster  and  Xlm  ^m'hm(^Ar+ii‘) 
as  the  held  scattered  by  the  cluster.  Of  course,  the  held  we  experimentally  measured  by  the 
procedure  outlined  in  Section  5.2  is  (r)  and  not  '^^Cm^rn{kN+ir). 

To  hnd  the  unknown  expansion  coefficients  of  the  experimentally  measurable  held,  we 
require  the  total  held  and  its  normal  derivative  to  be  continuous  across  the  r  =  a  interface, 
to  obtain 


®m'-fm(^iV+2®)  T  (-^m  T  -^m)  hfrn(^iV+2®)  —  (®m  T  Bm)  JmikN+lOj)  T  C rnH mik N +20j)  ^  (6.6a) 
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iJmi^N+20-)  +  (-^  m  +  A 

m  )  H'^{kN+2a) 

kN+i 


+  B^)  J'^{kN+ia)  +  CmH'^{kN+2a)  .  (6.6b) 


kN+2 

Eliminating  (B^  +  Bm)  from  (6.6),  and  using  (6.4)  with  (5.1a)  to  simplify  the  result,  yields 


Eliminating  Cm  from  (6.6)  and  then  using  (6.4),  (5.1a),  and 


(6.7a) 


{•Tm+l  (^Af+2®) )  kfm(kjV+l(^)J  {77m+l  (^A''+2®)  j  Jm(kjV+l(^)J 

{'7m+l  (^A''+2®) )  '7m(^7V+l®)  }  {77m+l  (^Af+2®) )  77m(/CjV_|_lO)  ]■  = 


to  simplify  the  result,  yields 


Bm  {hhn_|_l (/C7V+2®) )  7/m(/CjV+l®) } 


(6.7b) 


The  expansion  coefficients,  (B^  +  Bm),  of  the  wave  incident  onto  the  cluster  of  core  cylinders 
must  be  linearly  related  to  the  expansion  coefficients.  Cm,  of  the  wave  scattered  by  the 
cluster  because  the  Maxwell  equations,  material  properties  of  the  media,  and  the  continuity 
conditions  are  all  linear.  Traditionally,  this  relationship  is  denoted  by  the  so-called  T-matrix 
(transition  matrix)  [6-9].  Thus,  from  the  dehnition  of  the  T-matrix,  we  have 


Cm  =  Y.^. 


mn  V®n  +  B 

n) 


Combining  (6.7)  with  (6.8)  yields 


^  ^  i^^mn  Fmn)  -^n  Cm, 


where 


F  = 

J-  m.r). 


2i 


%Tmn  \_H n-\-\{k  N +2^) ,  77^  (/c^y+l®)  } 


Cm  — 


(6.8) 

(6.9a) 

(6.9b) 

(6.9c) 
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and  5mn  denotes  the  Kronecker  delta  function.  Let  us  compare  the  system  of  equations  for 
the  cluster  of  core  cylinders  (6.9)  to  the  system  of  equations  for  one  eccentrically  embedded 
core  cylinder  (5.12).  In  (6.9),  the  physical  and  morphological  properties  of  the  cluster  of 
core  cylinders  are  solely  contained  in  the  T-matrix.  Likewise,  in  (5.12)  the  effects  of  one 
eccentrically  embedded  core  cylinder  on  the  scattered  held  are  solely  contained  in  the  paren¬ 
theses  term^^  in  (5.12).  In  other  words,  the  T-matrix  in  (6.9)  plays  the  same  role  as  the 
parentheses  term  in  (5.12),  and  thus,  we  conclude  that  (6.9)  is  identicah^  to  (5.12).  In  light 
of  this  conclusion,  it  follows  from  Chapter  5  that  (6.9)  may  be  solved  via  the  Neumann  series, 
where  the  hrst  term  in  the  series  will  give  the  SCV  approximation.  The  SCV  approxima¬ 
tion  (an  orders-of-scattering  approximation)  and  its  physical  interpretation  were  presented 
in  Section  5.6  and  we  will  not  repeated  here. 

We  now  turn  our  attention  to  the  cluster  T-matrix  that  appears  in  (6.9).  If  we  denote 
the  T-matrix  for  the  ith  core  cylinder  as  then,  as  shown  in  the  next  section,  the  cluster 
T-matrix  will  be  composed  of  Th=^’  'N)  explicit  form  of  T^)  may  be  deduced  by 

inspection  of  (6.4a)  to  be 


{Jn+likN+lbi) )  Jmikibi)}  ^ 
(/uJV-|_l5j) ,  Jmikibi)} 


(6.10) 


6.2  Cluster  T-matrix 


To  derive  an  expression  for  the  cluster  T-matrix,  we  analyze  the  helds  inside  the  imaginary 
cylinder  of  radius  rmax,  which  contains  the  core  cylinders,  see  Figure  6.1.  The  scattered  held 
produced  by  the  ith  core  cylinder  is 

X]C'mV(/i^V+lPi),  \pi\>K  (6.11) 


^^The  Tmp  and  Tpn  quantities  in  yJ^pTmp^pTpnj  are  not  related  to  the  T-matrix.  This  rather  unfortunate 
notation  is  a  mere  coincidence. 

^^The  appearance  of  l/D^  in  (6.9b)  and  (6.9c)  is  not  important,  since  we  could  have  absorbed  it  into  the 
definition  of  via  Am  — >  IhmAm  as  was  done  in  Chapter  5. 


and  the  “external”  incident  field  on  the  ith  cylinder  is 

'^{Bra  +  Bm)^m{kN+ir).  (6.12) 

m 

According  to  Foldy’s  principle  [10],  the  effective  incident  field  on  the  Ah  core  cylinder  consists 
of  the  external  incident  field  and  the  scattered  field  produced  by  all  other  core  cylinders.  In 
other  words,  the  effective  incident  field  on  the  Ah  core  cylinder  is  given  by 

N 

C^'>'^{kN+iPj)-  (6.13a) 

m  j=l  m 

jA* 

Using  (6.2a),  with  w  =  r,  u  =  rj,  v  =  p^,  and  (6.2b)  with  u  =  r*  —  r^,  v  =  pj,  w  =  p^,  to 
express  \['m(/i^Ar+ii‘)  and  "^{kN+iPj)  in  (6.13a)  in  the  local  coordinate  system  of  the  Ah  core 
cylinder  yields 


E 


N 

B(i)  +  +  EE 

j=l  m 


expansion  coefficients  of  the  effective  incident  field 


$n(/CjV+lPi), 


where 


i  =  1, . . . ,  iV,  (6.13b) 


+  Bf  =  ^  (B^  +  Bm)  ^m-n{kN+lV^).  (6.14) 

m 

The  expansion  coefficients  of  the  scattered  field  produced  by  the  Ah  core  cylinder  are  related 
to  the  expansion  coefficients  of  the  effective  incident  field  via  the  T-matrix  for  the  Ah  core 
cylinder  and  thus,  we  have 


cf'  =  E^< 


(b 

ip 


N 


+ Bf  +  E  E 


j=l  m 


(6.15) 


where 


=  '^m-p{kN+i{rj  -  r*)). 


(6.16) 
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If  we  let 


B(2)  +  b(2) 

B(^-i)  + 

BW  + 

c(i) 

C(2) 


q(N-1) 


’xd) 

0 

0 

0 

0 

rp(2) 

0 

0 

0 

0 

0 

0 

0 

0 

.  .  .  rj^(N-l) 

0 

0 

0 

0 

rp(V) 

0 

0(1.2) 

O(dl^-l) 

QihN) 

0(2.1) 

0 

0(2.3) 

0(2, N) 

OiN-1,1) 

q{N-1,2) 

0 

q(N-1,N) 

OiN,i) 

OiN,2) 

OiN,N-l) 

0 

then,  in  the  above  block  matrix  notation,  (6.15)  reads 


(6.17a) 


(6.17b) 


(6.17c) 


(6.17d) 


[I-TO]e  =  TlB,  (6.18) 

where  I  is  the  identity  matrix.  Solving  (6.18)  for  C  and  using  [I  —  T0]~^  T  =  [T“^  (I  —  T0)]~^ 
to  simplify  the  result,  yields 


where 


e  = 


^  -  0 


-1 


(6.19a) 


(6.19b) 


We  will  hnish  the  derivation  of  the  cluster  T-matrix  shortly,  but  hrst,  a  few  computational 
remarks. 
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(i)  From  (6.10),  we  see  that  the  T-matrix  for  the  ith  core  cylinder  is  diagonal,  and  thus, 
so  is  7.  In  other  words,  the  computation  of  in  (6.19b)  is  trivial. 

(ii)  The  core  cylinders  must  be  labeled  in  such  a  way  so  that  the  dominant  (numerically 
large)  block  matrices  appear  near  the  main  diagonal  of  the  0  matrix.  This  is  im¬ 
portant  for  numerical  stability,  as  diagonally  dominant  matrices  tend  to  be  numerically 
stable  during  inversion. 

(iii)  It  is  tempting  to  compute  ^  by  expanding  [T~^  ~  0]  ^  in  the  Neumann  series,  i.e., 

[j-i  _  0]  =  T  +  TOT  +  TOTOT  +  ■  •  •  ,  (6.20) 

but  this  should  be  avoided.  The  right-hand  side  of  (6.20)  may  be  interpreted  as  an 
orders-of-scattering  approximation  and  thus,  a  large  number  of  terms  will  need  to  be 
retained  during  the  onset  of  localization.  Loosely  speaking,  localization  occurs  when  the 
wave  becomes  trapped  (undergoes  a  large  number  of  bounces)  inside  the  cluster  of  core 
cylinders.  However,  it  is  perfectly  sensible  to  use  an  orders-of-scattering  approximation 
outside  the  cluster  of  core  cylinders,  as  we  have  done  in  Section  6.1. 

Returning  to  the  derivation  of  the  cluster  T-matrix,  if  we  partition  the  ^  matrix  into 
block  matrices  and  denote  these  block  matrices  by  where  i  =  1, . . . ,  iV  and  j  = 

1, . . . ,  iV,  then  (6.19)  may  be  written  as 

N 

c®  =  ^  .  (6.21) 

i=i 

Substituting  (6.14)  into  (6.21)  yields 

N 

c«  =  ^  (B  +  B) ,  (6.22) 

1=1 

where 

Oim(ri)  =  ^m-n{kN+lYj).  (6.23) 
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The  expansion  coefficients  of  the  scattered  held  produced  by  the  ith  core  cylinder,  are 
related  to  the  C  coefficients  appearing  in  the  second  term  on  the  right  hand  side  of  (6.5d). 
To  establish  this  relationship,  we  express  (6.11)  in  the  global  coordinate  system  by  using 
(6.2c)  with  u  =  Fj,  V  =  Pj,  and  w  =  r  to  obtain 

^m-n{-kN+iri)C^^  Tn(A;Ar+ir),  |r|  >  |ri|,  i  =  (6.24) 

n  m 

Summing  (6.24)  over  i  yields  the  desired  relationship,  namely, 

N 

c  =  y^O«(-r,)C«.  (6.25) 

i=l 

Finally,  multiplying  (6.22)  by  r^),  then  summing  the  result  over  i,  and  using  (6.25) 

yields  the  cluster  T-matrix,  namely, 

N  N 

T  =  y^y^O«(-r*)^(*’^)0(^')(rj).  (6.26) 

i=i  j=i 

6.3  Connection  with  localization  and  random  matrices 

If  we  let  the  radius  of  the  host  cylinder  go  to  inhnity,  and  the  radius  of  each  core  cylinder 
shrink  to  zero,  then  the  scattering  geometry  shown  in  Figure  6.1  may  be  treated  by  Foldy’s 
method  [10].  In  this  extreme  limit,  we  can  treat  each  core  cylinder  as  a  dipole,  with  the  held 
scattered  by  each  core  cylinder  given  by 

k%^^PiG{r,ri),  i  =  (6.27) 

where  p*  is  the  dipole  moment  of  the  ith  core  cylinder,  and  G{r,  r*)  is  the  free-space  Green’s 
function.  The  free-space  Green’s  function  is  given  by  [11,  §2.2] 

G(r,ri)  =  ^i7o(Wik  -  g|),  (6.28) 

where  we  used  the  “negative  Dirac  delta  function”  convention  in  the  dehnition  of  Green’s 
function,  i.e., 

(V^  +  G(r,  r)  =  -(5(r  -  r). 
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To  find  the  unknown  dipole  moments  in  (6.27),  we  proceed  as  in  Section  6.2  by  considering 
the  effective  incident  held,  on  the  ith  dipole.  The  effective  held  incident  on  the  ith 

dipole  consists  of  the  external  incident  held  as  well  as  the  helds  produced  by  all  other 

dipoles;  thus,  we  have 

N 

r,),  *  =  1, . . . ,  iV.  (6.29) 

i=i 

j¥=i 

If  we  assume,  as  Foldy  did,  that  the  dipole  moment  is  proportional  to  the  ehective  incident 
held^^,  i.e.. 


V^  = 


(6.30) 


then  (6.29)  becomes 


N 


t;(efr)(r^)  ^  ^(inc)(r,)  +  V^),  ^  =  1,  .  .  .  ,  iV. 


(6.31) 


j=i 


In  principle,  the  proportionality  constant  Qi  is  arbitrary.  However,  if  we  assume  that  each 
core  cylinder,  as  well  as  the  host  medium,  is  non- absorbent,  then  conservation  of  energy 
mandates  that  Qi  must  satisfy  [12]  [4,  §8.3.1] 


=  exp(2iQ;j)  —  1, 


(6.32) 


where  ctj  is  the  phase  angle  of  g^.  Substituting  (6.28),  (6.32)  into  (6.31),  and  assuming  that 
Qi  is  the  same  for  each  dipole  (i.e.,  identical  core  cylinders)  yields 

FU^ff)  =  (6.33) 

where 


2ia  _  1 

F  =  I  G,  = 

2 

■uFff)(ri)' 

U(®®)(r2) 

,  = 

■t/(mc)(r^)- 

[/(mc)(r2) 

U("®)(r^) 

[/(ini 

(6.34a) 


^^This  assumption  forces  us  to  treat  each  core  cylinder  as  an  isotropic  scatterer.  It  should  be  noted  that  it 
is  possible  to  extend  Foldy’s  method  to  anisotropic  scatterers,  as  was  shown  by  Martin  [4,  §8.3.3]. 
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and 


Ho{kN+i\ri  -  Tjl)  for  i  j 

0  for  i  =  j 


(6.34b) 


From  (6.33)  and  (6.30),  we  see  that  if  F  has  a  zero  eigenvalue,  then  is  non-zero  even  if 
■jj(mc)  _  Q_  Qiea2;ly^  state  corresponds  to  perfect  localization;  we  dehne  perfect  localiza¬ 
tion  by  vanishing  of  the  total  (scattered  -|-  incident)  time-averaged  energy  density  sufficiently 
far  away  from  the  cluster  of  core  cylinders.  It  has  been  proven  by  Rusek  et  al.  [12]  that  this 
perfect  localized  state  cannot  exist  for  a  hnite  N.  Although  it  is  impossible  to  have  perfect 
localized  states,  we  may  have  quasi-localized  states  characterized  by  almost  zero  eigenvalues 
of  F.  Through  extensive  numerical  studies,  Rusek  et  al.  [12-14]  was  able  to  demonstrate  that 
the  density  of  eigenvalues  of  G  cluster  near  negative  one  for  sufficiently  large  dipole  densities. 
In  other  words,  if  we  denote  the  eigenvalues  of  G  by  g,  then  Iie{g)  ~  —1  for  almost  any 
random  distribution  of  dipoles  with  sufficiently  large  density.  Of  course,  the  eigenvalues  of 
G  are  related  to  the  eigenvalues  of  F  by 


/ 


(6.35) 


where  f  denotes  the  eigenvalues  of  F,  and  from  (6.35)  we  see  that  a  may  always  be  chosen 
such  that  /  ~  0  for  R.e{g)  ~  —1.  The  choice  for  a  is  dictated  by  the  imaginary  part  of  g.  For 
example,  if  Im(5f)  =  0  or  —  1,  then  we  may  choose  Q;  =  7r/2or7r/4,  respectively,  to  obtain 
/  ~  0  (quasi-localized  state).  Furthermore,  Rusek  et  al.  numerically  demonstrated  that  if  the 
number  of  dipoles  is  increased  while  keeping  their  density  constant,  the  eigenvalue  density  of 
G  approaches  the  Re{g)  =  —1  boundary  fastest  along  the  Im(5f)  =  0  line  (resonating  dipole 
line).  It  is  interesting  to  note  that  very  little  is  known  about  the  distribution  of  eigenvalues 
of  non-Hermitian  random  matrices  such  as  G  [15].  In  fact,  papers  devoted  to  this  study  only 
recently  started  to  appear  in  the  literature  [16-19]. 

Returning  to  the  original  cluster  of  core  cylinders  with  hnite  radii,  we  see  that  (6.19) 
is  analogous  to  (6.33).  More  specihcally,  the  role  of  F  in  (6.33)  is  analogous  to  the  role 
of  =7-^-0  in  (6.33);  recall  that  T  ^  is  a  diagonal  matrix  that  depends  only  on 
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the  properties  of  the  core  cylinders  and  the  0  matrix  depends  only  on  the  relative  position 
of  the  core  cylinders.  To  the  best  of  our  knowledge,  there  is  no  formula  connecting  the 
eigenvalues  of  —  0  to  the  eigenvalues  of  and  0.  However,  the  discussion  in  the 
previous  paragraph  suggests  that  quasi-localization  is  most  likely  to  occur  when  the  core 
cylinders  are  near  resonance.  Therefore,  with  a  hxed  “resonant”  T“^,  we  can  study  whether 
there  is,  indeed,  a  signature  of  quasi-localization  in  the  spectrum  of  0.  It  is  important  to 
stress  that  our  model  allows  us  to  do  this  numerically  and  experimentally,  unlike  the  above 
mentioned  Foldy-based  model. 

6.4  Summary 

In  this  chapter,  we  have  generalized  the  SCV  approximation  to  the  case  where  the  host 
cylinder  contains  N  eccentrically  embedded  core  cylinders.  This  was  accomplished  via  the 
cluster  T-matrix  approach.  We  demonstrated  that  the  SCV  approximation  retains  its  physi¬ 
cal  interpretation,  as  was  previously  discussed  in  Section  5.6.  Furthermore,  we  have  discussed 
how  our  model  can  be  used  to  study  localization,  and  suggested  some  avenues  for  future  re¬ 
search. 
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Overview:  Multiple  scattering  effects 

«  Develop  experimentally  verifiable  models 

►  Theoretically  sound 

►  Computationally  tractable 

►  AB  Millimetre  millimeter/sub-millimeter  Vector  Network  Analyzer 
9  Localization 


170  Teflon  &  Quartz  Layers  @  124.3  GHz 


X  (mm) 
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Overview:  Localization 
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What  can  we  measure? 


Alex  J.  YufFa  ayufFa@gmail.com  (MPL)  Modern  Electromagnetic  Scattering 


August  20,  2013  5  /  23  I 


Graf's  addition  theorems 

«  $n(fcr)  = 

►  ^m{kr)  =  y^$n(A:p)Qnm(ro) 

n 


►  '^mjkp)  =  y^^n(fcr)Qnm(-ro) 

n 

®  Onm(ro)  —  n(^ro) 

«  Regular  wavefunction: 

®  Outgoing  wavefunction: 

9  Mode-mixing 


r  >  To 
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Scattering  by  a  small  core  cylinder  inside  of  a  host  cylinder 


host  cyl. 


U^‘^\y\p)  —  +V^‘^\y\p) 

host  cyl. 

[/(3)(p)  =  V^^\p) 
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Host  cylinder 

U('"^)(r)  =  5^DA(r 

m 

u(sca)(j.)  ^  ^ 
m 

m 

9  Continuity  conditions  on  r  =  a 

|0“>  =  3 

{'Ati+I  j  Jm  (/C2a)} 


=U(2)  and  =  |:U(2) 


{Hm+i{kia)]Jm{k2a)} 

^ _ ] 

""  7vkia{Hm+iikia)] 

Jm  (/C2a)} 


•  Curly  bracket  notation 

V 

{fm+l{u);  gm{v)}  =  fm+liu)gmiv) - fm{u)gm+liv) 

u 
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Composite  cylinder 

m 

V^^\r-,p)  =  J2  {Bm^m{k2r)  +  Cm'^m{k2p) 

m 

V^^Hp)  ^J2^^^rn{k3p) 

m 

9  Continuity  conditions  on  r  —  a&L  on  p  —  b 
►  Graf's  addition  theorems 

^  ^  i^mn  Fmri)  -^n  —  Gm 


P  — 

i.J-  m.r).  — 


Tvkia^ 
2i  ^ 


^  ^  Omp{  Tq) Ap0p7j(ro)  I  {Hn+l{kltt)',  Hn{k2Ci}y 


iGm  —  ®m  I  ^  ^  Omp{  Tq)  Ap0p7i(ro) 


p 
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Composite  cylinder  cont'd 

9  Core  cylinder’s  contribution  is  small  contained  in 


{'^p+i(^3^)j  Jp{k2b^} 
{jp+i{k3h)-Hp{k2h)y 
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Physical  interpretation  of  the  SCV  approximation 

«  SCV  (orders-of-scattering)  approximation:  ^  —  ^ 

■^m  —  Gjji  —  ^  ^  I’o)^pOpn(l’o)®n 

n,p 


9  Physical  interpretation 

►  Opn(ro)®n  is  the  "screening”  effect  of  the  host  cylinder  on 

►  To)  is  the  “screening”  effect  of  the  host  cylinder  on 


Relative  error  in  energy 


0 _ 3 _ 6  9 


.... 

^0  = 

=  0° 

do- 

=  90° 

00- 

=  If 

>0° 

0  0 - ' — ^ ^ ^ 

0.2  0.4  0.6  0.8  1.0 


kslh 
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Host  cylinder:  Morphological-dependent  resonances 


•  Host  cylinder  ^  wavelength 

►  Ray  theory 


9  Caustic  radius 

►  Ray’s  angular  momentum:  \k2,0\rh 

►  mth  eigenmode  angular  momentum:  \m\h 
^  k2,rij'  —  r’caustic)  —  0 


m 


9  Bound  on  m 

►  Ray’s  angular  momentum:  |A:2|asin7 

►  Total  internal  reflection:  \ft\jV2  <  siny  <  1 

\ki\a  <  m  <  \k2\a 

9  SCV  approximation:  tq  <  rcaustic 

►  r’caustic  >  |A:i/A:2|  a  «  0.7a 
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Host  cylinder  resonance:  Energy  density 


9  Eigenmode 

►  m  =  228 

►  99.82385859...  GHz 


9  Caustic  radius 


►  r-caustic 


90° 


270° 
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Host  cylinder  resonance:  Spectral  radius 
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Scattering  by  a  cluster  of  eccentrically  embedded  cores 


host  cyl. 

host  cyl.  fmax<’'<« 
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N  core  cluster 

m 

y(iV+l)(r)  ^  ^  \Bm^^{kN+lv)  +  Cm-^^{kN+lv) 
m 

9  Physical  interpretation 

►  Cluster’s  incident  field 

U(^+i)(r)  + 

m 

►  Cluster’s  scattered  field 

m  (^7V+ir) 


9  Cluster  T-matrix 

►  Incident  and  scattered  coefficients  must  be  linearly  related 


Cm 


Tmn  (®n  +  Bn) 
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N  core  cluster  cont’d 

9  Continuity  conditions  on  r  =  a 
9  Definition  of  the  cluster  T-matrix 


mn 


An^G 


m 


^mGm  —  ^mTmnMn 

n 

^rnTrnn  {_Hn-\-l{kN-\-2(A)\ Hn{kN-\-l(A)'\ 

9  Same  linear  system  as  before 


C)mp(  i’o)^pC)p77,(ro)^  ^  Tj7 


►  Tmn  is  unknown 
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Cluster  T-matrix:  ith  core  cylinder 

o  Scattered  field  by  the  ith  core  cylinder 

m 

«  Effective  incident  field 


scattered  by  all  other  core  cyl. 

“external”  incident  field  / - -s 

- ^ ,  N 


E(»™  +  Sm)$„.(%+ir)  +  EE 

m 

9  T-matrix:  fth  core  cylinder 


j' 


j=l  m 


cf  =  Eri;’ 

N 

,  i  =  1, . . .  ,iV 

p 

j=l  m 

y(f)  _  _  {Jn+ljkN+lk)-,  Jmjkibi)}  (k^r^^(r■-r■)) 

{Hn+likN+M,JmikA)y^  r.jj 
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Cluster  T-matrix:  Block  matrices 


r  c(i)  1 

r  +  1 

=sr 

, - - s 

C(2) 

B(2)  +  b(2) 

e  =  [j-i 

e  = 

q{N-1) 

,  ®  = 

B(^-i) 

B(fV)  +  bW 

■t(i)  0 

0  0 

-1 

0  T(2)  0  0 

0  0  0  0 

0  0  rj.{N-l)  Q 

0  0  0  T(^) 


0 

0(1’2) 

OihN-1) 

OikN) 

0(2.1) 

0 

0(2.3) 

0(2.^) 

o(^-i.i) 

0(A^-1,2) 

0 

O(A^-l.N) 

o(^.i) 

0(^.2) 

0{N,N-l) 

0 
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Cluster  T-matrix:  Block  matrices  cont’d 


9  Global  coordinate  system:  pi  to  r 

N  N 

i=l  j=l 

►  block  matrices 


-J-(l.l) 

5^(1.2)  . 

^(2,1) 

5^(2, 2)  . 

.  ^(2,iV) 

^W2)  . 

.  ^{N,N) 

9  Compare  T  to  one  eccentrically  embedded  core  cylinder 
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Random  matrices  connection 


9  Foldy’s  method:  Isotropic  point  scatterers 

►  Assume  dipole  moment:  pi  = 

►  Conservation  of  energy:  =  2exp{2iai)  —  2 


=F 


UOnc); 


Gij  — 


Ho{kN+i\ri  -rjl)  if  f  j 
0  if  i  =  j 


9  Quasi-localized  state:  Eigenvalues  of  F  0 

►  Real  part  of  eigenvalues  of  G  — 1 

►  Choose  a:  Resonating  dipole 


Rusek  et  al.  PRE  51,  1995 
Rusek  et  al.  PRE  56,  1997 
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Outlook  and  Summary 
Outlook 

o  Role  of  F:  -  0 

Non-Hermitian 

T-i  diagonal  but  [T-i,0]  ^  [0,T-i] 

«  Role  of  G:  0 

«  Eigenvalues  of  random  matrices 

Localization  signature  in  the  spectrum 

Summary 

«  Produced  3  +  1  publications 
9  Created  start-of-th e-art  open-source  E(§zM  software 
9  Analyzed  interplay  between  absorption  and  causality 
9  Developed  experimentally  verifiable  model  to  study  localization 

Thank  you! 

Questions/Comments? 
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